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A. M. RODRIGUEZ,* P. A. LAGERSTROM,? ann E. W. GRAHAM? 


Douglas Aircraft Company, Inc. 


(I) ABsTRAct 


The paper deals with the problem of determining the geometry 


camber, twist, and mean angle of attack) which leads to mini 
mum drag of a wing of fixed plan form and given lift in super 
sonic flow. The linearized theory of an inviscid fluid is used 
A systematic study and extension is made of the method of 
orthogonal loadings. Several general theorems are proved for 
arbitrary plan forms and for plan forms with special symmetries 
These theorems will facilitate the practical determination of wing 
In the 


mathematical theory of Hilbert 


geometries of low drag Appendix, the relation of the 


methods used to the space is 


pointed out 


(II) SyMBOoLS 


local lift distribution 

local angle of attack 

Cartesian coordinates on base plane 

projection of lifting surface on base plane 

lift 

drag 

D/L? 

vector free-stream velocity 

subscripts referring to forward and reverse flow, 
respectively 

= notations for two types of reverse flow 

= coefficients in series expansion 


= average angle of attack 


(III) INTRODUCTION 


ii HAS LONG BEEN KNOWN! that the elliptical span 


wise loading gives the minimal drag for a three 

dimensional wing with fixed plan form and given lift 

in incompressible flow. The corresponding problem of 

finding the optimal loading and geometry in supersonic 
Received July 27, 1953 
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flow is considerably more difficult. In the latter case 
only a portion of the drag appears as kinetic energy in 
the Trefftz-plane. As a consequence, the Trefftz-plane 
methods of incompressible flow cannot be used, and it 
is necessary to take the chordwise loading into account, 
in addition to the spanwise loading. 

Recently, Jones* has established general criteria for 
optimal loading in supersonic flow and has constructed 
the optimal loading for special plan forms. However, 
his criteria are, in general, nonconstructive. Graham* 
and Sears* have independently developed methods for 
decreasing the drag of a given wing by combining vari 
Refer 


ence 3 shows the advantage of using loadings that are 


ous load distributions in suitable proportions. 
orthogonal-—that is, having zero interference drag. It 
is shown how orthogonal loadings may be constructed 
and how, by the use of an increasing number of orthog 
onal loadings, one may approach arbitrarily close to 
the optimal loading. 

This report contains a systematic study and exten 
the 
theorems regarding orthogonal loadings are proved. 


sion of method of reference 3. Several general 
In particular, it is shown how the criterion for muini- 
mum drag of Jones is a direct consequence of these 
theorems. Some theorems are also proved for special 
plan forms with certain symmetries. 

The theorems proved here are valid, in general, only 
if the leading-edge suction is neglected. The signifi 
cance of this is discussed in Section (VIT). 

In the Appendix it is pointed out that the methods 
used here actually are a special case of the general 
methods of the theory of Hilbert space. The remainder 
of this paper is, however, independent of this Appendix 

In principle, the methods of reference 3 may be used 
to determine a loading that is arbitrarily close to the 


optimal loading. For a given plan form, one may start 


co! 


U 
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with, say, the flat-plate loading. An arbitrary second 
loading may be orthogonalized and combined with the 
flat-plate loading so as to decrease the drag for a given 
lift. The process may then be repeated with a third 
loading, etc. Examples of this method were given in 
In practice the computational difficulties 
The theo- 


reference 5. 
become prohibitive after the first few steps. 
rems given here will facilitate the choice of orthogonal 
However, further research is needed to 
A method for estimating how 


loadings. 
shorten the procedure. 
much larger the drag of a given loading is when com 
pared with the minimum drag would be desirable. 
This is an important problem for further research. 
Other important problems were mentioned in reference 


) 
» 


Dr. R. Sedney has assisted the authors in the prepa- 


ration of this report. 


(IV) ORTHOGONAL LOADINGS 


We shall consider three-dimensional wings of zero 
thickness. 
then of equal magnitude and opposite sign at corre- 
sponding points of the upper and lower surfaces. For 
convenience we shall denote, by p(x, y), the difference 


The pressure distribution on such a wing is 


in pressure between the lower and the upper surface 
i.e., the local lift distribution. The local angle of at 
tack will be denoted by a(x, y). 
local downwash divided by the free-stream velocity. 


It is equal to the 


For a given plan form, a(x, vy) determines the twist and 
We shall 
refer to these features as the geometry of the wing. 
Note that the “geometry” in this sense is changed by a 
By the formulas of linearized wing 


camber and also the mean angle of attack. 


rigid rotation. 
theory, p(x, vy) determines a(x, vy) and, conversely, 
a(x, vy) determines p(x, y). 

A lift distribution p and its associated angle of attack 
distribution @ will be referred to as a load distribution 
or a loading (p, a). By linearity, it follows that, if ¢ 
is a constant and (p, a) is a load distribution, then 
(cp, ca) is another load distribution. Two load dis- 
tributions are said to be of the same type if one is a 
constant times the other. For convenience we shall 
write 

SS f(x, v) dx dy = rr LS 


wing 
With this notation 


lift =L= fpdsS (la) 
drag = D= f[padS (1b) 
Note that, as mentioned in the introduction, any lead- 
ing-edge suction is neglected in the expression for the 
drag. We define 


d = D/L? (2) 


and it follows from Eq. (1) that d depends only on the 
type of loading. 
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For a wing fixed in space, we shall denote the flow 
with vector free-stream velocity U as forward flow and 
the flow with vector free-stream velocity — U as reverse 
flow. The drag is the component of the aerodynamic 
force in the direction of the appropriate free-stream 
velocity, and the angle of attack is measured with 
respect to the appropriate free-stream velocity in each 
case. Witha given loading (/, a) in forward flow will 
be associated two loadings in reverse flow, (p, &) and 
(p, &). In the first case the lift distribution is kept the 


same 


Pp = p, & determined from p (3) 


In the second case the geometry of the wing is kept the 


same. Then, 


a = —a, fp determined from 4 (4) 


Frequent use will be made of the following reciprocity 
theorem:’ If (f,, ay) is an arbitrary loading in forward 
flow and (p,, a,) is an arbitrary loading in reverse flow, 


then 
S pra. dS = Sp a, dS (5) 


In particular, if (p;, a1) and (fo, ae) are two loadings in 
forward flow, then Eq. (5) may be applied to either (/,, 


a) and (ps, &) or (pi, ay) and (ps, de): 


S pr& dS = S poo dS (6a) 
Spree a5 = S pow ds (6b) 


By linearity, two loadings (fi, a1) and (ps, a2), with 
drags D, and V2, respectively, may be added to form a 
new loading (f1 + f2, a1 + a2). The drag of the new 


loading is 
S (pi + pr)(a, + ay) dS = Di + D2 + Dy 
where Dy» is called the interference drag. The inter 
ference drag is the same if the flow is reversed accord 
ing to Eqs. (3) or (4), 
S (pra + pom) dS = J (prom, + pride) dS = 
S (Poe ot Pde) ds (7) 
The first equality is proved as follows: 
S (pre + pea) dS = S (Pree + pra,)dS, by Eq. (3) 
= S (pom + pia)dS, by Eq. (6a) 
= S (pom + pi&)dS, by Eq. (3) 


The second equality is proved in the same manner. 


Following reference 3, we define (/;, ai) and (ps, as) 
to be orthogonal if 
Dy = S (pres + Pea dS = 0 (S) 


It follows directly from Eq. (7) that, in this case, 


S (pres aa ped) dS a S (Pree + Pod) ds = 0 (9) 


Furthermore, if ); = ps = p, Eq. (7) implies that 
D=D=D (10a) 
L=L (10b) 


st 
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In case @ constant, 


L L 


but, in general, LZ is different from 1 and —L. 
A set of orthogonal loadings, 


(Pry A) = (Pr, a1), (Po, ae), ... (11) 


is called complete if an arbitrary loading (p, a) may be 
expressed as a linear combination ef the elements of this 


set -1.e., 


p= > Appr, a= } B Ay, (12) 
k=] k=l 
There are infinitely many such complete sets. One 
may start with an arbitrary collection of loadings, 
orthogonalize them (cf. reference 3), and add additional 
loadings to make the set complete. Given such a set 
of loadings, we define 


L. = SpdS, De = S pro, dS (13) 


Then, by linearity and orthogonality, from Eqs. (la 
and (1b), 


> tele, D = > a,2D, (14 
] 


k k=] 
We shall also use the notation 


d, = D,/L,? (15) 


Expansion (12) is analogous to a Fourier series (cf. 
Appendix), and the coefficients a, may be found in a 
similar manner as the coeflicients of a Fourier series 
namely : 

Theorem 1.—With the above notation, 


7 


I ; 
a, = 2D, / (pa, + pra) ds 
This is easily verified by substituting the expansions 


(12) into the integral and using orthogonality. 


Theorem 2.—With the above notation and assump 
tions, (Py, &) and (/,;, a) are two complete sets of load- 
ings in reverse flow. The expansion of (p, &) and (f, 
a), respectively, in terms of these sets have the same 
coefficients a, as (p, a) has in terms of (p;, a,). This 
follows easily from Eq. (7) and Theorem 1. Note that, 
with an obvious notation, Eq. (10) implies 


L, = L,, 


D, = D, = D, (16) 
(V) Minimum DRAG 


The loading for a given plan form which gives the 
This 


property depends only on the type of loading. A 


minimum drag for a given lift is called optimal. 


necessary and sufficient condition that a type of loading 
be optimal is thatd = D/L*isa minimum. As pointed 
out above, the quantity d depends only on the type of 
loading. We shall now show that, if a complete set 
of loadings is known, the optimal loading may be con- 


nono snr 
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structed. The optimal loading is denoted by (/,, 


It has an expansion in terms of the complete 


Qopt,)- 


set (11 


Daas Ss a,.p , Ss 
Las “KE KY f , a 


l 4 l 


dy, 


where the a, are to be found from the condition that d 
log D 2 log L, the 


necessary condition for a minimutn is 


is a minimum. Since log d 


(1/D) (OD /oa, (2'L) (OL /Oa 0) 


Hence, from Eq. (14), 
(1/D)(2a,D,) (2/L)(L 0 
or 
a, = (D/L) (L,/D 


which can also be written as 
a [L/>> (1/d,)| (L,/D 


It can be shown that these values of a, actually vield 
a nunimum, and, since we obtain a unique solution for 
the a,, the optimal load distribution is unique. There 
fore we have proved: 

Theorem 3.—In terms of a given complete set of 
loadings [(p;, a,)] the optimal loading is given by 


Popt (D/L) > (L/D) p (17 
l 

Qop. = (D/L) & (L,/D,)a Is 

L = (D/L 2d (Li*/D (19) 

D = (D/L)? >> (L,2/D (20) 

Ld = D/L = L/)>, (1/d 21 


Note that only the lifting members (L, # 0) of the set 
occur in the above formulas. This is necessary because, 
if one adds a zero lift loading (whose local lift is not 
everywhere zero) to the set, the total lift remains the 
same while the drag increases. Furthermore, all the 
contributions of the individual loadings to the lift have 
the same sign as L. That this is necessary is clear be 
cause, if for an arbitrary set of loadings L is positive 
and one adds a term a,L, which is less than zero, the 
drag increases while the lift decreases. To bring the 
lift back to the value L, one must add a positive term 
that makes the drag still larger. 

For the reverse flow, with the complete set |(p);, @ 
corresponding to Eq. (11), the same type of calcula 
tions used in proving Theorem + can be used to prove 

Theorem 4.—The loading (pp, Popt.» Bopt,) 18 opti 
mal in reverse flow, and its expansion in terms of 
[(De, &)] has the same coefficients as the corresponding 


expansion in forward flow. That is, 
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Fic. 1. Characteristics of the optimal loading. 














(Deoos.): = (Dent,); (29) 

Popt, = (D/L) > (L/D) Pr (233) 
k=1 

Bop, = (D/L) D> (Le/Dr)& (24) 
k=) 


where the subscripts f and r refer to forward and re- 
verse flow, respectively. |For the case in which the 
complete set [(p,, &)] is used, see below. } 

For an arbitrary p, from the complete set, consider 


the following expression : 
(1/2D,) S” Pe(copt, + Hops.) dS (25) 


From Eq. (6a) and Theorem 1, Eq. (25) can be written 


as 
l ; | 
oD. (Pr Cpt. + Popt, &) aS = 
(1/2D, )S (pe Qopt. + Pope. Ok) dS 
= A, 
= (D L) CE, D,) 
= (D/L) (1/Dy) Sp, aS 
Therefore, 


S Pp, [Qopt, i Qopt. a (2D bias == {() 


From Eq. (12) it follows that the result is 
If the quantity @ = [a@pr. + 


for all px. 
true for any pressure p. 
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ope (2D L)]| is regarded as an angle of attack dis 
tribution in forward flow and /p is the corresponding 


pressure, then 
S padS = 0 


Since the drag is zero only for the zero loading, it follows 
thatw = 0. Thatis, 

Qopt, + Bp = 2D/L (26 
at all points of the lifting surface. This is the criterion 
given by Jones* for the optimal loading of a lifting sur 
face having total lift L. 

Formulas similar to Eqs. (23) and (24) cannot be 
proved in general when the set [(/,, @)] is used. In 
general, L, # 1,. 

The following formulas may be proved for general 
plan forms: 

Theorem 5.—The necessary and sufficient condition 
that a loading (p, a) be orthogonal to a flat-plate load 
ing (Pi, a, = constant = c) is that L = L. 


Proof: 


i = cS p aS = S ap ds 
cL = S ap is = S ep; dS = — S ap, ds 


Hence, 
c(L — L) = S (ap + ap,) dS 


which proves the theorem. 
For the flat plate itself, LZ) = —/,.° 
an orthogonal set [(~;, a;)] where the flat-plate loading 


Consider now 


is the first member a; = constant = c. For the reverse 


flow, form the set [(p;*, a;*)] which is identical with 
(pi, a) ], except that a;* = —a = c so that p,* = 
— py. Since (pi*, ay* 


tiplicative constant (—1), [(p,;*, a:*)] is a complete 


) and (/;, a) differ only by a mul- 


orthogonal set by Theorem 2. Furthermore, by 


Theorem 5 and Eq. (16) 


L,* = L,, D,* D, = D, (27) 


II 


Consider now the plan form in reverse flow with the 
same lift L. The drag for optimal loading must then 
be the same as in forward flow, D. The coefficients of 
the expansion of the optimal loading in terms of [(p,*, 


a;*)] are 
(D/L) (2,°/D,") = (D/2£) {L./D) 


Hence we have an alternate way of forming the optimal 
loading for reverse flow. Comparing this result with 
Eqs. (17), (1S), (23), and (24) and using the uniqueness 
of the optimal loading, we have: 

Theorem 6.—Using the notation of Theorems (3) and 
(4) the optimal loading in reverse flow may be expressed 


as 
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Popt.r = (P de ache = 2 Orbe = QI aebr* 
, emg 
(Qopt.), > > aa. = >» (Ly. & 
f 1 P 1 


but, in general, 
(Qupt )ys x ( Qopt 


Let us now consider the relation between (a,,,;), and 


We have 


Qop: )y More closely. 


(Qont Jp = QA)Q) + S Q;,Q; 
7 au 
as * S * 
(Qopt Jy = AQ + a, 
k=2 
= Qa; _ QA); 


Thus the optimal loading in forward flow may be de 
composed into a flat-plate distribution a, = dja; = ay 
and an orthogonal distribution 


> 


azn = % Ay) 


The shape corresponding to ag has the property that it 


has the same lift in forward and reverse flow. Then, 
(Qopt.)s = Qa + ap (28 
(Qopt)y = AA — Ap (29) 


where the a, and a z in reverse flow are measured with 
respect to the reverse flow. For the corresponding 


pressure distribution, 


(Pp, mdf = (Popt as 
= pa + Pr 
= — ps = Pp 


These conclusions are illustrated graphically in Fig. 1. 
Note that a, is constant for all sections. 

Conversely, if one has a flat-plate distribution a, 
and an orthogonal distribution ag such that the asso- 
ciated pressure distributions satisfy 


pa + pe = — px + Pp 


then the loading (a4 + ag, Ps + par) is optimal. This 
is so since, if the same pressure distribution is used in the 
reverse flow, then the shape in reverse flow is ay — ap. 
Hence the sum of the angle of attack distribution ‘is 


2a, = constant in accordance with Eq. (26). 


(VI) SpecraAL PLAN ForMS 


Some theorems applicable to plan forms of particular 
shapes will now be proved: 

Theorem 7.—For plan forms with all edges super- 
sonic, no side edges, and a straight trailing edge, a 
loading (p, a) must have zero lift if it is orthogonal to 


WINGS OF FIXED PLAN FORM 





Y 
STRAIGHT 
SUPERSONIC 
U TRAILING EDGE 
or oe __y 


SUPERSONIC —— 
LEADING EDGES 


Angle of attack distributions that offer no drag reduc- 
tion 


Fic. 2 











Fic. 3. A plan form symmetrical about the y-axis 


the flat plate (p;, a; = c) and if a is constant along lines 
parallel to the trailing edge. 
Proof: If the shape is unaltered and the flow is re 


versed the loading is two-dimensional —.e., 
pb = Ka = —Ka 
Pp} = Kea = —Kay = Ke 
By the reciprocity Theorem (5) 
L= fpdS = -(1/0c)S pe dS 


= —(] oS pra a5 = Kf a dS 
L= Sp dS = KS a ds ~K fa dS 
Hence L = —L. it fol 
lows that L = 0, which proves the theorem. 
From the remarks following Theorem 3 it follows that, 


Since, by Theorem 5, L Ls 


for the plan forms considered in Theorem 7, the flat 
plate drag cannot be improved by adding the type of 
angle of attack distributions as described in Fig. 2. 
Next consider a plan form of the type shown in Fig. 3 
Theorem 8.—For a plan form symmetrical about the 
y-axis and a load distribution |p(x, y 
(a) If a(x, vy) = a(—x, y) and the load distribution 


» aly, V } 


is orthogonal to the flat plate, then its lift is zero. 
(b) If a(x, vy) = —a(—x, y), then the load distribu 
tion is orthogonal to the flat plate. 


(c) The optimal geometry has zero twist. 


Proof: From the symmetry of the plan form, it fol- 
lows that in case (a) L = —L. But since (p, a) is 
orthogonal to the flat plate, L = L by Theorem 5. 
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Hence L = 0. Incase (b), L = L from the symmetry 
of the plan form and the antisymmetry of the angle of 
attack distribution. Orthogonality to the flat plate 
then follows from Theorem 5. 

Part (a) of the theorem shows that, for the plan forms 
considered, the drag of the flat plate cannot be decreased 
by combining it with orthogonal loadings whose angle 
of attack distribution is symmetrical about the y-axis. 

Consider now an optimal shape a(x, y) for the sym- 
metrical plan form of Fig. 3. By symmetry, the angie 
of attack distribution in reverse flow a(—.x, y) must be 
optimal. Applying this to the decomposition (28), 
(Qopt ry = Gas + ap, it follows that ay(—x, y) + ap 
(—x, y) must be optimal in reverse flow. However by 
Eq. (29) the optimal shape in reverse flow is a4(x, y) — 
ap(x, y). Since a4(x, y) = ay(—x, y) = c, it follows 
that ag(—x, y) = —ap(x, y). 
wing corresponding to such an ag distribution is sym- 
Thus, for a plan 


The geometry of the 


metric with respect to the y-axis. 
form symmetric about the y-axis, to improve the flat- 
plate loading it is necessary to use antisymmetrical 
angle of attack distributions (symmetrical geometry). 
For the optimal geometry the mean angle of attack at 
each spanwise station is that given by the flat plate 
component, which proves case (c). 

For computational purposes, the following equations 
are convenient: 

If the plan form is symmetric as above and aq and a» 
are odd in x, then the interference drag between the 


loadings (a), Pi) and (ae, ps») is 
Dy = 2S apo = 2S ap, (30) 


This follows from the symmetry and the reciprocity 
relations. 

Furthermore, for a symmetric plan form, let ZL; and 
D, be the lift and drag, respectively, of the flat-plate 


distribution a; = |. Let a» be odd in x and let a = 
€,a; + Ca be the optimal combination of a; and a2 with 
total lift L equal to L;. (This is not necessarily the 
optimal a for the plan form.) Then, 

Did D 

7 Lan ® 


Hence, if a is the average value of a over the plan form, 
D/D; = & (31) 


Thus the relative drag reduction is simply 1 — a. 


(VII) REMARKS ON LEADING-EDGE SUCTION 


Most of the above theorems are valid only if the 
leading-edge suction is neglected. Orthogonality may 
be defined even for the case of leading-edge suction. In 
this case the drag should be defined not by Eq. (1b) but 
by a formula obtained by applying the momentum 
theorem to a region bounded by a control surface en- 


The interference drag would then 


closing the wing. 
be defined as the drag of a; + as minus (drag of a; + 
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(b) FINITE PRESSURE 
ACTING ON FINITE 
AREA 


(c) LOWER PRESSURE 
ACTING ON LARGER 
AREA 


The approximation of leading edge thrust with non- 
singular loadings 


Fic. 4 


drag of as). Two loadings would be called orthogonal 
if the interference drag is zero. 


sion in orthogonal loadings could then be defined as 


The method of expan- 
before. However, the theorems based on the reciproc- 
ity relation would no longer be valid. In particular, 
the drag of a wing is, in general, altered by reversal of 
the flow direction if there is leading-edge suction. 

The theorems of the present paper may then be inter- 
preted as follows: They are obviously valid for the 
case of wings with supersonic leading and trailing edges 
and with side edges parallel to the flow. If some 
edges are subsonic, they would be valid if it is assumed 
that leading edge suction is not realized. However, 
there is also the following possibility. Assume that 
the optimal loading for a given plan form is determined 
taking leading-edge suction into account. There are 
indications that the geometry of the optimal loading 
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inay require such a warping of the subsonic leading edges 
that leading-edge suction does not occur even theo- 
retically. f 

Consider a plan form having subsonic leading edges 
and a twist and camber distribution such that leading 
edge singularities are produced. It should be possible 
to approximate as closely as desired the lift and drag 
of this loading, using a modified twist and camber 
chosen to eliminate the singularity. 

The modification would be made by warping the ex- 
treme forward portion of the airfoil so that the flow 
crosses the leading edge smoothly. Something equiva- 
lent to leading-edge thrust is then produced by a large 
but finite negative pressure acting on a small but finite 
area that is inclined forward. (See Fig. 4.) 

If the preceding point of view is adopted, it then 
becomes unnecessary, from a practical point of view, 
to study loadings that have leading-edge singularities. 
Either the class of loadings considered may be re- 
stricted, or (as in this report) we may assume that 
leading edge thrust does not exist. In the latter case, 
loadings with singularities become highly undesirable, 
since by warping the extreme forward portion of the 
airfoil to eliminate the singularity something equivalent 
to leading-edge thrust can be obtained. The loadings 
with singularities, being undesirable, would ultimately 
be eliminated in the drag reduction process and so are 


of no concern 1n minimum drag problems. 


(VIII) APPENDIX 


It will be indicated below how the methods used in 
this paper are related to certain standard methods in 
the theory of Hilbert space, as described in references 
Hor 7. 

The projection of the wing plan form on the x-y-plane 
is a compact two-dimensional point set with a Euclidean 
measure. Consider the measurable functions on this 
set. These functions will be denoted by a(x, y) and 
interpreted intuitively as angle of attack distributions. 
The corresponding pressure distribution, p(x, vy), 1s 
a certain linear integral operator, 
The Hil- 


bert space will consist of those functions a(x, y) for 


found by applying 
given by linearized wing theory, to a(x, y 


which the corresponding p(x, v) exist and whose drag 
is finite. Functions that differ only on a set of measure 





+ Some of these ideas were suggested by R. T. Jones. 
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zero will be identified. The inner product between two 
angle of attack distributions, a; and a», will be defined 


as 
[a, Qe | = (1 2) JS (Piae + Poa ds 


which is half the interference drag D,. between the two 
distributions. It is easily verified that [a), a2] is bi 
linear and symmetric. The norm of a function, a@ , is 


then given by 


al\? = fa, al S pa dS D 


} 


VD 


It is zero only when a@ vanishes 


Since the drag is always positive, the norm || a 


is real and positive. 


outside a set of measure zero. One may also prove 
Schwarz’'s inequality, 
[ay, a] S iia a 


The function space with inner product defined above 
can then be seen to satisfy all the axioms of a separable 
Hilbert space. The mappings a > & and a > a de 
fined in Eqs. (3) and (4) are then algebraic isomor 
phisms of the space, and from Eq. (7) it is seen that the 
mappings are also metrically isomorphic. The concept 
of orthogonality defined by Eq. (S) is now identical 
with the concept of orthogonality from the point of 
view of Hilbert space. The expansion of an arbitrary 
function in terms of orthogonal loadings as given by 
Eq. (12) is then the ordinary expansion of an element in 


Hilbert space in terms of orthogonal base vectors 
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SUMMARY 


The transient distribution and the thermal 


stresses in an idealized wing structure considered by Hoff and 


temperature 
Torda in reference 1 are determined. Only the effects of aero 
dynamic heating and of heat conduction are included; radiation 
and convection effects are neglected. The present work differs 
from that of reference 1 in that the conduction from the cap to 
the web is considered when the temperature of the cap is calcu 
lated, and the spar cap temperature is assumed to be a function 
of both space and time. Graphs of temperature and thermal 
stress distributions are presented, and the results are compared 


with those of reference | 


List OF SYMBOLS 


) 


a, b = equivalent lengths for the web and flange, Fig. 2 
1, B = constants in the Laplace transformations in Eqs 
(17) and (18) 
Pas = area of the web 

A, = area of one flange 

A in = Fourier coefficient for 7) in Eqs. (84) and (86 

B.. = Fourier coefficient for 7) in Eqs. (33) and (35 

ri = thickness of cover plate, Fig. 1 

Cc = specific heat; constant; contour 

C,(t) = function of time for web, Eq. (44) 

C,(t) = function of time for flange, Eq. (44 

d = thickness of web, Fig. 1 

E = Young’s modulus 

h = coefficient of heat transfer from boundary layer to 
wall 

H = one-half of height of web, Fig. 1 

K(t) function of time, Eq. (47) 

L = width of cover plate, Fig. | 

M = number of roots of Eqs. (28) and (29) 

p = complex variable in Laplace transformation, Eqs 
(8) and (9) 

r? = p/k = —a? 

q? = (h + p)/k = +f? 

t = time 

To = initial temperature in cap and web 

7; = web temperature 

T> = cap temperature 

T, = boundary-layer temperature 

YN = Laplace transform of 7(x, ¢), Eq. (8) 

W = d, Fig. 1 

We = 2c, Fig. 1 
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distance measure from center of web, Fig. 2 


coefficient of thermal expansion 


a —_ 

Gm, Bm = roots of transcendental Eqs. (25), (26), (28), and 
(29 

K = thermal diffusivity 

a = thermal stress in web 

g, = thermal stress in flange 

p = density 


v and fas subscripts denote partial differentiation with respect 


to these variables 


INTRODUCTION 


ee AND EXPERIMENTS agree in predicting high 
temperatures for the boundary layer surrounding 
the wing of an airplane flying at supersonic speeds. 
The full effect of these temperatures upon the behavior 
of the wing structure is not yet known. It is the pur 
pose of the present paper to solve some of the problems 
involved through a rigorous calculation of the distribu 
tion of the temperatures and the thermal stresses in a 
model of a supersonic wing. <A typical cross section of 
a wing is shown in Fig. 1. 

It is assumed that the upper and lower cover plates 
of the wing are connected by a number of shear webs. 
The analysis presented here deals with a representa- 
tive portion of the structure consisting of one shear 
web and one-half of the length of the coverplates be- 
tween adjacent shear webs. It is assumed that the 
boundary-layer temperature on the upper and lower 
surfaces has the same given constant value. Because 
of the temperature differences between the boundary 
layer and the cover plates, the latter are heated. At 
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Fic. 1. Cross section of the skin-shear web configuration. 
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the same time heat is conducted away from the cover 
plates by the web. Asa result of this heat conduction, 
the temperature of the cover plate near the juncture 
of the web and cap is considerably lower than the tem- 
As the heat 
flow from the boundary layer into the coverplate de 


perature halfway between the two webs. 


pends upon the temperature difference, the determina 
tion of the time-space history of the temperature is 
rather complex. 

The problem has been solved by the method of the 
Laplace transformation. For any fixed time after the 
beginning of heating, the calculations yielded tem- 
perature distribution curves that are very similar to 
curves obtained from experiment.’ The stresses calcu 
lated are equal to, or greater than, the yield point of 
the material after heating (or the equivalent super- 
sonic flight) of SO sec. 

The results of this work are in satisfactory agreement 
with the solution of the same problem by N. J. Hoff 
and Paul Torda in reference 1, in which much more 
far reaching simplifying assumptions were made. 


FORMULATION OF THE PROBLEM 


Hoff and Torda, in Appendix I of reference 1, treated 
the problem of the aerodynamic heating of an idealized 
spar shown in Fig. 1, subject to certain assumptions 
listed on pages 104-105 of reference 1. These assump- 
tions will be retained here. The most important as- 
sumptions state that (a) aerodynamic heating is caused 
solely by the heating of the spar cap through the adja- 
cent boundary layer, (b) the effects of radiation and 
convection are ignored, and (c) the wall temperatures 
and the heat-transfer coefficients are the same at the 
top and bottom wing surfaces. 

The aerodynamic heating in (a) is proportional to 
the difference in temperature between the boundary 
temperature and the spar cap temperature 7%(x; 1) 
In addition to the assumptions listed by Hoff and Torda, 
their treatment of the problem was simplified by the 
following assumption: The temperature of the spar 
cap was assumed to be a function of time only. This 
temperature was determined on the assumption that 
the heat loss by conduction from the cap to the web 
was negligible. In this way the simple analytical re- 
sult that the cap temperature 7» is an exponential func- 
tion of the time and is independent of any space co- 
ordinate was obtained. To determine the web tempera- 
ture 7)\(x; ¢), the usual equation for the one-dimen- 
sional conduction of heat was solved under the assump- 
tion that the end of the web was kept at the known tem- 
perature 7.(f). Since 72(¢) is a simple function of the 
time ¢, the solution for 7\(x; ¢) can be determined 
easily (reference 1, pages 106-108). Thus this as- 
sumption leads to an approximate solution of the prob- 
lem. 

The problem will be treated as one involving a ther- 
mal coupling of two bars rather than as two simplified 


single bar problems. The cap temperature 7>2(x; ¢) 
is now a function of both space and time, and the heat 
loss to the web by conduction will not be neglected in the 
calculation of the cap temperatures. 

For convenience in the mathematical formulation and 
solution, the boundary-layer temperature 7) will be 
Thus 7); < 0, 7> 


receives heat from the boundary layer and if the initial 


assumed to be zero. 0 if the cap 


temperature 7) is less than 7,. The differential equa- 
tions are 


(73)¢ = «(7))-, (GO <zs<s@4 (1) 


(73); = «(T2)er — AT: (€ Sx SO (2) 


The two bars are assumed to be composed of the same 
material. The additional term (—/7%) in Eq. (2) rep- 
resents the rate of heat flow into the cap from the 
boundary layer and is positive since 7, < 0. A more 
detailed treatment of the heat balance is given in refer- 
ence 4. 


The boundary conditions are 


(7), =0 v= 0 (3) 
T; = 7; (x =a (4) 
w(7}); = w(ls), (x=a (5) 
(T2), =0 (x =) (6) 


Eqs. (3) and (6) stipulate that no heat flow takes place 
across the ends of the composite bar. Eqs. (4) and 


(5) represent the continuity conditions at the junction 
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Fic. 2. Simplified representation of Fig. 1. (The points are 


located at the quarter points of the web and flange, respectively. 
The values of a and + refer to the cross section for which results 
are shown in Figs. 4, 5, 6, 7, and 8 
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7 the web and the cap; w/w» is the ratio of the cross- 
sectional areas of the two bars (Fig. 2). 
[he initial conditions are 


7, = Ty = constant| 
T. = Ty) = constant} 


[he mathematical problem to be treated is specified 
by the differential Eqs. (1) and (2); the boundary con- 
ditions of Eqs. (3), (4), (5), and (6); and the initial 
conditions of Eqs. (7). The solution will be obtained 
by the I_aplace transformation method.’ 


SOLUTION BY LAPLACE TRANSFORMATION 


Ir 7 = 7\(x; t), then the Laplace transform of 7 will 


he denoted by 7*, where 


T* = T*(x: p) = | e T(x: 2) dt (dS) 
0 


The Inversion Theorem (reference 2, pages 244-245) 


then states that 
Piss ft) = (1/298) / eT *(x; p) dp (9) 
« Y i 


over a suitable contour in the complex p-plane. 
Che application of Eq. (8) to Eqs. (1) and (2) yields 
the ordinary differential equations for the transforms 


of JT, and 7%, 


d?7,*/dx*) — r*T;* -To/ (10) 
d*T.* dx?) — q°T2* -7To k (11) 
where 
y- = p/K, q- (h + Pp). k (rz) 
have been introduced for convenience. The trans 
formed boundary conditions are 
d7,*/dx = 0 (x = 0 (13) 
i.” = 7.” (= a (14 
w,(71* ws 1 9* (yr ae @ (15) 
al." /dx = @ (x = b (16) 


Solutions of Eqs. (10) and (11) which satisfy Eqs. (13) 
ind (16), respectively, are 

1,;* = A cosh (rx) + (T0/p (17) 

TT,” = B cosn [g(x — b)]) + [To (2 + p)] C8) 


To satisfy the two remaining equations, Eqs. (14) and 


15), it is necessary that 


A cosh (ra) + (7) p) = B cosh [q(b — a)] + 
[70 (A + p)] (19) 


7A sinh (ra) = —wqB sinh [g(b — a)] (20) 


The solutions of Eqs. (19) and (20) are 
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A = (To P) (ht Kh) (Wo Ww) (1 q) x 
sinh q(b a){l A(p 4 | 


B [7 (h + p)] (hex) [11 /r) sinh (ra)] [1 Alp 


where 


A(p) = r sinh (ra) cosh |q(b a)| + 


W2/w)qg sinh |g(b a)]| cosh (ra 23 


The quantities y and g are functions of p, by Eq. (12 

The final values of the transforms 7\* and 72* are 
given by Eqs. (17) and (18), with the values of A and B 
as given by Eys. (21) and (22). To evaluate 7, and 
T, by Eq. (9), the expressions for 7)\* and 7 :* must be 
examined in detail 

The transforms 7)\* and 7>* are single-valued ana 
lytic functions of the complex variable p, since the ap 
propriate power series expansions in Eqs. (21), (22), 
and (23) contain only integral powers of p. The ap 
propriate contour C is thus the usual one (reference 2, 
page 244) shown in Fig. 3 

The integral from y to to y +120 in Eq. (9) de 
termines the temperature. If the integral over the 
circular portion of the contour C approaches zero as the 
radius tends to infinity, then the evaluation of Eq. (9 
is equivalent to the evaluation of the residues of the 
integrand. The proof that the integral over the curved 
portion of C tends to zero follows the procedure usual 
with such problems (reference 2, pages 367-369) and 
need not be given here. 

The poles of Eqs. (17) and (18) must be considered. 
First, consider the poles of Eqs. (17) and (18) which are 
not the result of the vanishing of A(p). In Eq. (17) p = 
—h (that is, g = 0) is not a pole, and in Eq. (18) p = 0 
(that is, r = 0) is not a pole, since in each case the inte 
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grand is finite. However, p = 0 is a pole of Eq. (17) 


and p = —hisa pole of Eq. (18). In each case a simple 

calculation shows that the residue is zero. The only 

residues remaining are those corresponding to the 

zeros of A(p). The proof that all zeros of A(p) are 

real and negative can be carried through in the usual 

way (reference 2, pages 266-267) and will be omitted. 
If p > O (that is, r > 0, g > 0), then A(p) > 0 and no 

If p < 


positive roots exist. 0, infinitely many roots 


exist. Two cases must be considered. 
Case 1. Let p < O and fA + p < O-—that is, assume 
p < —h. This assumption states that r? and g? are 


both negative in Eq. (12). Set 
—-e =r = p/x, —B = q’ = (A+ p)/wx (24) 


rhe vanishing of A(p) is then expressed by the vanish 


ing of the transcendental equation, 


a sin (aa) cos [B(b — a) | + 
%2/w:)8 sin [B(b a)| cos (aa) = Q (25) 
where 
a® — Bp? = h/k (26) 


Eqs. (25) and (26) have infinitely many real roots, all 
of which are simple; a@ and 6 are both real (a? > h/x). 
rhe roots can be located approximately by slide rule, 
and the results can be refined by Newton's method. 
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a? + B? = h/k 29 


For purposes of numerical calculation, Eq. (28) can be 


rewritten in the form 
a tan (aad) = (w./w,)8 tanh |oi/ a 30) 


<. h/n, ( 5 


30 is 


From Eq. (29) it follows that 0 < a? 
8? < h/«. 
As ft increases in magnitude, the number of 


Thus the number of roots of Eq. 
finite. 
such roots steadily increases. In a practical case the 
roots can be found by the same procedure as used in 
Case 1. 

The discussion of Cases | and 2 assumed that # > 0. 
If h < 0, which would correspond to the case of heat 
loss to the boundary layer, then the second case above 
(p <0, h + p > O) is impossible, and these additional 
roots would not appear in the solution. Thus the twe 
cases h > Oandh < Oare quite different mathematically 
only the case h > 0 is of interest in the present work 

The roots of Eq. (25) determine an infinite set of 
terms, and only cosine terms will arise from these roots 
However, the finite number of roots determined by Eq 
30) will yield cosine terms for the web and hyperbolic 
cosine terms for the cap. The entire set of functions 
vields the complete set for the problem 

To evaluate the residues, (d/dp){A(p)| must be evalu 
ated subject to the condition that A(p 0. It is 


convenient to simplify the expressions for 7) and 


Case ?.Let p < 0 but A + p > O—that is, assume 
-h <p <0. This assumption states that r? is nega before evaluating these residues explicitly. After the 
tive and g’ is positive in Eq. (12). Set necessary manipulations and simplifications have been 
| __ made, the results can be written in the form 
—aQ =r = p/x, B° = @ n-+p)/« (27) 
fio : : hans oe > . 
Ihe vanishing of A(p) is then expressed by the vanish- I(x; t) = 2, Am COS (a,X)E 
. ° F ” ey 
ing of the transcendental equation, 
ul 
@ sin (aa) cosh [B(b a)| — T(x; t) = ) B,, cosh 8,(b — aje “° ~ 
uz %,) B sinh [B(b — a)]| cos (aa) = 0 (28) é 
2 B,, cos 8, (6 — aje “* * 32 
where Ws 
In Eqs. (31) and (32) a@,, and 8,, (m = 1, 2,..., 17) are the roots of Eq. (30), which are finite in number. For 


these terms, Eq. (32) contains hyperbolic cosine terms. 
are the roots of Eq. (25), which are infinite in number. 


plifications, are 


The quantities a,, and 8,,(m = \J + 1, WM + 2, , 
The equations defining the Fourier coefficients, after sim- 


Bi = 
oT 
5. jcosh 8, (6 — a) COS @,,A ' ; 
cosh 8,,(6 — a) + : (a) (b — a)(we wi) | m = 1, 2, _M 
(h/x) sin a,a ' COS a,,A cosh 6,,(6 — a) 
4.3 
A,, = B,,|cosh 8,,(6 — a)/cos a,,a] 34) 
ek 
B.. = witli 
5. * jcos B,,(b — a) COS Q@,,f j 
cos 8,,(6 — a) — : (a) + (b — a) (w2/w { 
(h/x) sin a,,a | COS a@,,l cos B,,(6 — a) 
(m= M+1,M+2,..., rc 35 
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Fic. 4. Web and cap temperatures (7°) as a function of the time (¢) for various values of .Y 
, B,, [cos Bn(b — a) cos Amt | (36 
Eqs. (31) and (32) represent the solution; the coefficients A,, and B,, are defined by Eqs. (33) to (36). The con 


stant initial temperature is 7); the temperature of the boundary layer is zero. 


TR 


ANSIENT TEMPERATURE DISTRIBUTION 


The same numerical values used by Hoff and Torda 
in reference 1 will be taken. The numerical constants 


are: 


0.0186 sq.in. per sec 
60°F. (initial temperature ) 
600°F. (boundary-layer temperature 


é.0 


= ().011623 per sec. = a Cpe 


4.6875 in 

16.6875 in 

heat-transfer coefficient from the boundary layer to 
the wall = 90 B.t.u./(hour ft.? °F. ) 

specific heat of steel = 0.1406 B.t.u./(Ib. °F.) 

density of steel = 489.6 Ibs. per cu-ft 


= ().375in 


The transcendental Eqs. (30) and (25), respectively, 


may be written as 


Qm tan (4.6875an,) — 7.5(0.62493 — amn*)’* X 


tanh [12(0.62493 — a@,,2)"*] = 0 


n= 1,2, ey A) CZ) 


a, tan (4.6875a,,) + 7.5(a,,” — 0.62493) “? x 


tan [12(a,? — 0.62493)*] = 0 
(m= M+1,..., ©) (388) 


Eq. (37) has two roots (J = 
eae 


and the first twelve roots of Eq. 


finitely many roots (m = 


Eq. (37) (m = 


> 


(38) (m = 3,... 


mi Am 
l 0.32243 


2 0.78613 


3 0.823801 
} 0.92298 
5 1.01509 
6 1. 18050 
rf 1.31349 
8 1.51716 
9g 1.66522 
10 1.77901 
1] 1.99656 
12 2.22193 
13 2.35014 
14 2.49799 


The quantities a, 


a = A ” 


m 


Ey a 


7 A 


Eq. (38) has in- 


The two roots of 


, 14) are given in Table 1. 


0 
0) 
0 


Q).2 


TABLE | 

p Am 

72178 1.50170 
08325 — 1.42586 
22894 0.33038 
17640 0. 25644 
63677 0.33716 
SOS15 0.04046 
04897 0.00117 
29493 —().00772 
16562 —(). 05023 
59372 —(). 00946 
83339 —(Q OOO138 
07655 0.00331 
21319 0.01662 


36960 


_and b,, are defined as 


(—7»), 


b,, = B,, 


h 


0. 0000309 
0.79156 
0. 26989 
—0.11439 
0.07318 
—(). 02327 
0.00116 
0.00530 
—0. 00798 
0.00458 
—0.00013 
—(). 00186 
0.00238 
—() 0O101 


— 7%) (39) 


Eqs. (31) and (32) become, in this case, after suitable 


adjustment to the initial temperature of 600°F., 


T(x; t) = 600 


40 


m 


> Um 


l,2 


Kant 


COS (Q@p,x)e- 


(40) 
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for 


whi 


T\ 





36 


on 


; of 
iq. 


0) 


} 


THERMAL 


STRE 


SSES IN A WING 13 





- WwEB 


BOUNDARY ty LAYER TEMPERATURE - 


FLANGE i 


T 





















- t ol 
| — | : _— ! 


| 














[oo } { + } 
| isankell 
& 400-—-+ et Bt thn eS 
uj a 
w i } 
P | 
<q 300+ +—_——+} \y 4 
a 
- = — =) ~— —+—_—_——+ + + || 
_ 2S | | 
~ 200-— = \ = 
— | | 
: eS Se OS ee ee | 
EET 
100} } t | | | : 
INITIAL TEMPERATURE — 
aii rT vee | 
| | | | | | | | 
Oo t i 1 1 i | i Bas A, J 
O | 2 3 4 5 8 ‘ 17 
DISTANCE IN INCHES FROM CENTER OF WEB 
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. = , (40) < ) bec - equs ( ity < = (), Cc 
T2(x;t) = 600 — 540 j b coak Ole ~ bie" 4 4 ind (41 yecome equal to unity it / the con 
be vergence is poorest at ¢ = 0. If 14 terms are taken in 
> ' — Eqs. (40) and (41), the temperatures can be computed 
nes ?m COS Bm(X — Oe { (41) to the nearest tenth of a degree at 15 sec. Asi 
: increases, fewer terms are needed until ¢ > 400 sec., 
Eqs. (40) and (41) can now be used to plot the tem- after which only two or three terms are needed. In 


peratures 7\(x; /) and 72(x; t) as functions of x and of ¢. 
The curves of Fig. + show the variation of the tempera- 
ture as a function of the time 
The nine points selected correspond to the ends and the 
quarter points of each bar, Fig. 2; thus .; corresponds 
to the junction of the cap and the web. 
perature at the junction (x = a) may be computed from 
7\(a; t) or from 75(a; ¢), this calculation was used to 
check the arithmetical work. 
and xy were drawn as one curve, since the temperatures 
This curve is the 
same as the curve shown in Fig. 1, page 
91, which was computed on the basis of zero heat con- 


for various values of x. 


Since the tem- 


The curves for xs, X7, As, 


were the same to within 1°F. or less. 


7 of reference 


duction into the web. 
ror’ = 


be - 


for example) lead to the values at the end of the bar, 


T2(b; 


the results 


0.99536, 0.99944 (42) 


14 
>» .,, = 
m l 


T,(0; 0) = 62.5°F., 0) = 60.3°F. (43) 


which should be compared with the initial temperature 
Ty = 60°F. Since all exponential factors in Eqs. 


most heat conduction problems it is difficult to com 
near ¢ = 0. The usual 
alternative 
for small 


pute temperatures accurately 


procedure in such cases is to derive forms 
for the solutions which 
values of ¢ (reference 2, page 274); the results in the 


converge rapidly 
present problem would be too cumbersome to use. In- 
stead, it is simpler to compute additional terms until 
the required accuracy is obtained. As ¢ increases, 7; 
and 7, approach the boundary-layer cemapeteiate of 
600°F. 

The graphs in Fig. 
perature as a function of the distance from the center 


5 show the variation of the tem- 
of the web for various values of the time. 


THERMAL STRESSES 


The thermal stresses in a thin plate can be calcu- 
lated from the equation derived in reference 3, pages 


399-401: 
—akT + C 


c= 


where o is the normal stress perpendicular to the cross 


section of Fig. 1, a is the coefficient of thermal expan- 
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sion, and C is a constant that may depend on the time. 
[he actual bar temperature may be used for 7’, since 
any added constant will not change the stress; this re- 
sult will be apparent later. 


In the present problem the normal stresses in the 


web and in the flange may be denoted by 


eS = akT, oie Cat 


ee ee com) 


where the subscripts w and / refer to the web and to the 


va 
\| 


. 


where 7, = 
are given by the equations 


O,(V;t) = 


THE AERONAUTICAL SCIENCES 


= 7, and 7) = 7%, in agreement with the previous notations of Eqs. (40) and (41). 


—ak7T\(x; t) + akK(t), 


TANUARY, 1954 


flange, respectively. At the junction of the web and 
the flange, 7, = 7); since the stresses are equal at the 
same point, 
Cy = C, = akKtt (45 
where A(t) is a function of the time alone 
The vanishing of the stress resultant 1n the cross sec- 


tion of the spar is expressed by 


"a "bh 
d / Cy AX + 4 J g,ax =0 4H) 
« 1 Zia 


Eq. (46) vields the following value of A (¢ 


[1 + 2(A,/A,)] (1 a) / Ti (x; t) dx + (2A,;/Aw) [1/(6 a) | J Ts(x; t dx, 47 


The thermal stresses 


ox; t) = —akT72(v; t) + aEK(t) 48 


Phe two integrals that appear in Eq. (47) are simply the average temperatures of the web and of the flange, respec 


tively 


« 


*h » 


l. 


‘hese average values follow immediately from Eq. (40) and (41 


(1. a) / T\(x; t) dx = 600 - 340 | 


Lm ee. 


: the results are 


—_ 
> 


» An(1/a,a) sin (a,a)je “™ | (49 


1)] / T(x; t) dx = 600 — 540 > bd, [1 Bn(b — a)] sinh B,,(b — aje~**"" + 


a)| sin B,,(b ono 50 


The average temperatures, Eqs. (49) and (50), can easily be tabulated numerically, and from these tabulated values 


A(t) can be computed. 


The thermal stresses in the web are plotted in Figs. 


Hand7. To avoid overlapping graphs, Fig. 6 shows the 
results for ¢ = 50, 100, and 150 sec., and Fig. 7 shows 
the results for ¢ = 200, 300, 400, 500, and 600 sec. 


The thermal stresses in the flange are too small to be 
shown conveniently on these graphs, but some informa 
tion about these stresses can be obtained from Fig. 8, 
which shows the variation of the thermal stresses with 


time for various values of x. 


DISCUSSION OF RESULTS: COMPARISON WITH THE 
RESULTS OF HOFF AND TORDA 


The curves of Fig. 4 indicate the slowness with which 
the temperature of the web increases with time. After 
10 min. of heating, the center of the web is at 350°F. 
The junction of the web and flange (x;) heats much 
more rapidly, and the temperature reaches 585°F. in 
10 min. or only 15°F. below that of the boundary 
layer’ At points some distance from x;—namely, the 
pots Xs, X7, Xs, and xy—the temperatures are very 
The curve for these points coincides with 


nearly alike. 


The thermal stresses are then determined by Eqs. (48) 


the curve of reference |, Fig. 7, which was based upon 
the assumption that no heat flow by conduction takes 
place in the web. The assumption made by Hoff 
and Torda in reference | would make the temperature 
at x; equal to that on the curve (x6, 7, Xs, and x»); thus 
the end of the web is assumed to be at a higher tem- 
perature than is actually the case. This difference 
amounts approximately to 30°F. in the range from ¢ = 
100 to 300 sec. This difference is small because of the 
dimensions chosen for the cross section. The ratio 
2A, A, = 19.2 in this case indicates that the flange 
area is much larger than the web area. Thus less heat 
will be conducted away by the web in this case than if 
the ratio 2A,;/A, were smaller. If this ratio were 
smaller, the simplifying assumptions of reference | 
would lead to serious errors for the web temperatures. 


The curves of Fig. 5 show how pronounced is the 
effect of a large ratio of flange area to web area. The 
temperatures along the flange are very nearly constant. 
At ¢ = 200 sec., for example, the temperature at the 
end of the flange is 30°F. higher than that at the junc- 
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Fic. 7 


f 


reference | is equivalent to taking the ratio 2A ,/A,, as 


In effect, the assumption « 


tion of web and flange. 
infinity. This ratio is 19.2 in Fig. 2; from Figs. 9 to 
12 of reference | it is easy to see that the results of the 
present calculations should be in close agreement with 
this simplifying assumption. The curves are con 
tinuous everywhere but have discontinuous first de- 
rivatives at the junction of the web and the cap. The 
discontinuity in the first derivative follows from Eq. 


(5), since we/w, = 7.5. Thus 





rhermal stress (kips,/in.? 





in spar web (1 kip = 1,000 Ibs 


dT \(a; t)/dx = 7.5dT>(a; t)/dx 


for the curves shown in Fig. 5. 

The graphs of the thermal stresses shown in Figs. 6 
and 7 show the same trend as those shown in Figs. 9 
to 12 of reference 1. The stresses are always maxi- 
mum at the center of the web. The tensile stresses at 
the end of the web are quite small; the stresses in the 
flange are not shown because of their small magnitude 
The largest compressive stress in the flange is 3,500 Ibs 


per sq.in. 





16 JOURNAL OF THE AERONAUTICAL 
90- ,--y + Ty T T t 1 T T T T ] ] 
ee SS 
80: 1 


70; 


60- 


50 


40- 


30 


4. 
STRESS (KIPS/IN. ) 


20; 





+ 


os | ‘- 





05 “1 j00 | 200 ' 300 1 400 
| 4 t (SEC.) } | } 4 


+ + + 


Fic. 8. Thermal stresses (kips/in.*) as a function of time for 
various values of x (1 kip = 1,000 Ibs.). 


The variation of the thermal stresses with time is 
shown in Fig. 8S. The maximum stress occurs at the 
center of the web at ¢ = 225 sec.; the maximum value 
is 80,300 Ibs. per sq.in. The stress decreases slowly 
At the end of 10 min. the stress is still over 
A typi 


cal curve for the flange is also shown (xy); this curve 


with time. 
45,000 Ibs. per sq.in. at the center of the web. 


shows clearly how small the compressive stress is in the 


flange at all times. 


SCIENCES JANUARY, 1954 

The maximum stress computed in reference 1 occurs 
at the center of the web at t = 200 sec. (Fig. 10); the 
stress = 80,000 Ibs. per sq.in. This result is in excellent 
agreement with the present result. 


In conclusion, it may be stated that when the flange 
web ratio is sufficiently large, the temperature (and in 
particular the thermal stresses) may be computed with 
sufficient accuracy by the simplified treatment of refer- 
ence 1. If the flange-web ratio is small, the more ac 


curate treatment of the present paper should be used. 


The method used to solve this problem can be ex- 
tended directly to the case of asymmetrical heating 
that is, different boundary-layer temperatures at the 
top and bottom surfaces. In addition, the case of a 


flange of variable thickness can be considered. 
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+ The title of the graph of Fig. 10, page 94, of reference 1 is in 
error. The time ¢is 200 sec. and not 100 sec. as indicated. Simi 
larly, Fig. 11 of reference 1 should read ¢ = 100 sec. instead of 


t = 200 sec 
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A General Method for Solving 


Problems of 


the Unsteady Lifting Surface Theory 
in the Subsonic Range 


H. G. KUESSNER* 
Institute for Fluid Dynamics and Applied Mathematics, University of Maryland 


SUMMARY 


An oscillating lifting surface of any shape with any downwash 
The Kutta condition of finite ve 
locity at the trailing edge can be satisfied in a general way only 
The 
a member of the family of orthogonal 


distribution may be given 


if one introduces appropriate orthogonal coordinates 


lifting surface must be 
used by 


surfaces. Corresponding wave functions have to be 


which the wave equation can be separated. By means of these 
functions a characteristic function G may be constructed, the 
normal derivative of which at the lifting surface has the proper 
ties of a Dirac 6-function. The kernels of the integrals, which 
represent the exact velocity potential and the pressure, are com 
posed of this function G and its derivatives. Examples are given 
for the elliptic and circular plate and for the infinite long ribbon, 
The solution of the 


for which the wave functions are known 


two-dimensional case is also given 


(1) INTRODUCTION 
ie THE THEORY OF STABILITY AND OF FLUTTER, one 
needs aerodynamic derivatives for very small har- 


Meas- 


urements of these derivatives being difficult and not 


monic oscillations of the airplane or parts of it. 


suitable, one depends mainly on the linearized lifting 
surface theory in order to obtain them. These theo- 
retical derivatives yielded reasonably good values for 
the critical velocity of stability loss and for the damping 
in many cases compared with experimental values. 
Although boundary-layer effects have certain influences 
mainly for oscillating control surfaces, the trend of the 
functions obtained is represented very well by the 
lifting surface theory for nonviscous fluids, applying 
Kutta’s condition. 

In the supersonic range one can obtain the velocity 
potential from the downwash on the lifting surface in a 
simple way, using Kirchhoff's formula and the point 
solution of the wave equation, because all perturbations 
act only downstream. Therefore it has been possible 
in the last years to work out some aerodynamic deriv 
atives of the oscillating lifting surface for all wing shapes 
of interest, for instance, rectangular, triangular, and 
sweptback wings; but this work is not yet completed 

In the subsonic range, although every airplane must 
pass it and most of them will never surpass it, little has 


Received June 16, 1953 
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forschung, Goettingen, Germany 


been completed until now, except the two-dimensional 
case and the case of Mach Number zero. Therefore 
this subsonic range deserves a new treatment. 

If the pressure on the lifting surface is given, one ob 
tains without difficulty a general integral representation 
for the velocity potential, using the well-known point 
But if the 
and this is 


solution of the wave equation at rest. 
downwash on the lifting surface is given 
mostly the case—then the point solution is not appli 
cable to the corresponding integral representation of the 
solution because it is not possible to satisfy the Kutta 
condition with it in a simple manner. The experience 
on solving the two-dimensional problem has shown that 
this can be done in a general manner only by using ap- 
propriate orthogonal coordinates. The lifting surface 
has to be a member of the family of the orthogonal sur 
face system, and the wave equation has to be solved 
by products of appropriate wave functions. 

but there are only a few wing shapes, for which or 
thogonal coordinates are known—circular and elliptic. 
Less is known concerning the wave functions of these 
coordinates. Either we have to restrict ourselves to 
these shapes, or we have to use approximate methods, 
or we have to develop new orthogonal coordinate sys 
tems and new corresponding functions. Because ap 
proximate methods are in many cases cumbersome 
and uncertain, it seems correct to put the problem on 
the strong shoulders of the mathematician and to de- 
velop new wave functions. 
the two-dimensional case has been 
McLachlan! 


One needs correspond 


The solution of 
achieved mainly by using the book of 
on Mathieu's wave functions. 
ing books on Lamé’s functions for the circular plate, on 
for the 


the wave functions—not considered by Lamé 


elliptic plate, and on other wave functions. Brief re- 
views of the existing knowledge of these functions are 
given by Whittaker and Watson,” by M. T. O. Strutt,* 
and by E. W. Hobson.‘ As detailed knowledge of these 
functions is not available to the author at present, only 
an outline of the theory is given in this paper. Ex- 
amples containing the basic equations are given below 
for the elliptic and circular wing shapes and for the in- 
finite ribbon. 

Some time ago, the transonic range 
8 = 1—has been considered as a limit of the linearized 


Mach Number 
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lifting surface theory, giving infinity. But Heaslet, l l 

Lomax and Spreiter,® and also N. Rott® have shown X =Ix; Y= V1 — 8 ly; Z = V1 — 8 lz (2 
that, for finite reduced frequencies w of the oscillating , 

lifting surface, the derivatives of the two-dimensional o = lv: exp (ivf + wB*x)-p*(x, y, 3) 3 
case are finite and fit very well into the curve as func- 

tions of 8, as shown by I. E. Garrick.’ The infinity p — Po = pw: exp (tvt + wB*x)-p*(x, y, 2) (4 


at w = 0 may be a resonance effect, which vanishes also crs ; ; 
a : ‘ : where X, ¥, Z are common cartesian coordinates; / is 
for a finite span of the wing.® 


: i : a characteristic length, for instance, the half main chord 
For the stability theory one needs very low fre 


of the wing; 


: an g; v is the oscillatory frequency; v is the ve- 
quencies near zero... To overcome the above-men- : = ee ; 
phe ; : locity of the undisturbed flow in the direction of the 
tioned difficulty at 8 = | in another way, it has been 2s j : : 

: ; et ink, . positive x-axis; cy) is the velocity of sound, and py the 
suggested to start the linearization of the problem of i ‘ S ; : 
a : . ; density of the undisturbed fluid; p is the pressure; and 
small oscillations from the exact nonlinear differential er : ; : St 
: : : ape: Sige ¢ is the perturbation velocity potential. @ and p* are 
equation of the velocity potential for a finite, but small, rs rid ; 
ge : : Sanaa ’ dimensionless, as are w, x, *, y, and g. 
airfoil thickness. But Heaslet, Lomax and Spreiter,° 


and also Busemann’ have shown that for finite airfoil The substitutions of Eqs. (2) to (4) involve the so- 
thickness a steady potential solution as starting point called general Lorentz-transformation of coordinates.'° 
of the linearizing process does not exist at 8 = 1 unless With these substitutions, Euler's equations can be 
one admits surfaces of discontinuity (shocks). There- transformed into 
fore, for very small values of w near zero, we have to 
use the three-dimensional lifting surface theory. (. 4 Oo or x, y, 2) + p*(x, ¥, 2) = 0 (5 
Ohe 
(2) Tue Basic EQUATIONS 
(V2 + x«?)o*(x, y, 2) = 0 (6 
In order to obtain simple formulas, we use the fol- 
lowing substitutions for a harmonically oscillating lift- (V2 + «?)p*(x, y, 2) = O (7 


ing surface: = ; : he 
Ihe wave equations, Eqs. (6) and (7), are the same as 


8 e tvl vl (1) for resting sources. Therefore, we may use the known 
= : o= — K = ? ae ° 
; co(1 — B?) solutions of it. 


We can assume either the pressure p*(x, y, + 0) or the downwash 
w*(x, y, 0) = $:*(x, y, 0) 
on the lifting surface. Consequently, we have two kinds of solutions for the velocity potential in the exterior, 
o*(x, v,2) = SS K(x, x’, y, vy’, 2)p*(x’, y’, + 0) dx’ dy’ (8) 
o*(x, v,2) = SS K(x, x’, y, y’, s)w*(x’, y’, 0) dx’ dy’ (9) 


The integrals have to be taken over the lifting surface. If we insert these integrals into Eq. (5), we obtain the 


pressure 
* ss we , , , , ce) a , , 
p*(x, vy, 2) = —S S p*(x’, v’, +0) dv’ dy’ (a+ > RAS, x 3 ¥ 4S (10 
x 
* 7 ae as >. ; 
p*(x, v,2) = —S Sw*(x’, y’, 0) dx’ dy’ (wo + > Ase. VS (11 
x 


Both kernels A and K have the difference property 
K(x, x", y, y', 2,3’) = K(x — x’, y — y’,2 — 2’); (s’ = 0 (12 
Therefore, we can write all four solutions, Eqs. (8) to (11), in two different forms, Eqs. (Sa), (9a), (10a), and (11a), 
for instance, 
o*(—x’, —y’,s) = SS K(x — x’, y — y’, z)w*(—x, — y, 0) dx dy (9a 
The kernel A is the same in Eqs. (9) and (9a). But if we introduce other coordinates and transform the integrals, 


the integration process of the two solutions, Eqs. (9) and (9a), may be different. 


l is 
rd 
ve- 
the 
‘he 
nd 
ire 


as 


vn 


UNSTEADY LIFTING 


According to Eq. (12) we have 
| 


OK /Ox = —OK/On’ (13 


Cherefore, we may replace the operator |w + (0/0”) 
in Eqs. (10) and (11) by [w — (0 Ox’ 
calculate, especially if we introduce general orthogonal 


This is easier to 
coordinates into Eq. (11), since x’ holds on the lifting 
surface. 

Eq. (11) or (11a) is mostly needed for the calculation 
oft forces and moments of the oscillating lifting surface. 
But Eq. (9) or (9a) is also needed for the theory of sta 
bility of airplanes. The main lifting surface of the air 
plane induces a varying downwash on the tail surface, 
the induction of the tail surface to the main surface 


being negligible. 


In order to obtain convergent integrals for A,, if z = 0, 
we substitute, by means of the wave equation, Eq. (6), 


O° F 0° O° 
= —_— (x: + : — ) F (1S 
O02” Ou? Ow" 


following a suggestion of I. E. Garrick. 


(4) CONDITIONS FOR THE KERNEL A 
According to Eqs. (6) and (9), the kernel A has to 
satisfy the wave equation 
(V?+ 7)AK = 0 (19 
Furthermore, we must apply Eq. (5). Multiplying it 
with the operator 


[w + (0/Ox)] (0/dz (20 


we obtain 


R) 
(w As ) pete, v3) = 0 (21 
Ox i 


The general trend of the functions p*, ¢,*, and @,.* is 


represented in Fig. 1. Therefore, we have, on the lift 


ing surface (z = 0), 


Inserting Eq. (22) into Eq. (21) and integrating the 
linear difterential Eq. (2]) for z = 0, we obtain the 


SURFACE THEORY 19 


Eq. (8) has been used for iteration methods and ap 
proximations in order to invert Eq. (8) into Eq. (11 
without the knowledge of the kernel K 


3) THe KERNEL A 
The kernel K is well known 


solutions of the wave equation at rest in the subsonic 


Using the known point 


range, 
F(u, 2 Hy? (kV u? + 3? 14 
exp “Vue +wu-+s2 


F(u, w, 2 15 


and integrating the equation, Eq. (5), we obtain the 


kernels for the two- and three-dimensional case,! 


fe) > 
u)da Hy? (k Vai + 2? (16) 
Oz 
O exp (—t V a? + w? + 2? 
da 17 
Oz Vai t+ w?4+ 2? 


pressure gradient on the upper side of the lifting surface, 


. 


x 
(a lim p.*(x, y,z) = —w? lim $:*(x’, v, 3 
z=-+0 z + Of - 
exp w(x’ — x) dy’ 0 (23 
At x = —©, all perturbations must vanish, but not 
at « = +o, because of the wake. Let x(y) be the 


coordinate of the leading edge. From the Eqs. (9) 
and (23) we have 


(b lim / A(x — x’, y yd 
eS oe -Q (24 





“ y + 


above and below 
the lifting surface 


outside the lifting 
surface 


every where 


Fic. 1 Trend of the functions p*, ¢:*, @:* 








In order to reproduce the given downwash on the lifting 


surface, 


o.*(x, y, 0) = w*(x, y, 0) ( 
by Eq. (9), we have the necessary conditions 
; jO, if O or w 0 (26) 
(c) A,(u, w, 0) = me 
lo, ifu=w=0 
97 


“<é) 


(d) / / K(u, w, 0) dudw = | 
The function A,(u, w, 0) is a so-called Dirac 6 function 
or a zero function, also playing an important part in 
the theory of the Bessel and Mathieu functions. 

The Kutta condition prescribes a finite velocity at the 
trailing edge of the lifting surface. Therefore, we con 
struct the kernel 


l 


Laur 


= exp (—7xr) 


(g) lim K 


r= Vu? + w? + 2? (30 


This is the so-called retarded solu 
In the two- 


and its derivatives. 
tion, needed for obvious physical reasons. 


dimensional case, Eq. (30) goes over into 
lim AK = (2/2)Hy? (xr) | 

; (30a 
r= Vu? + 2’ | 


(5) ORTHOGONAL COORDINATES? 


In order to determine the kernel A according to the 


above-mentioned conditions, we must introduce an 


orthogonal coordinate system (u, v7, w) such that both 
sides of the lifting surface become a special surface of 


; ; ss the orthogonal surface system. Let 
(e) AK = A(G+ CG’) (28) 
; ; a x = X(uU,v,W); y = y(u,v, WwW); 2 = o(u,v,w) (31) 
by means of a regular solution G of Eq. (19), giving al T 
ways finite velocities, and by a singular solution G’ of With the auxiliary functions 
Eq. (19), giving infinite velocity only at the leading U-2? = (dx/du)? + (Oy/du)? + (d2/du)? 
edge of the lifting surface. We prescribe, for the de- d 
rivative, J = (Ox/Ov)? + (Oy/dv)* + (O2/d0v)* (32 - 
(hatte: =~ 9 9 = 4.2) W-? = (Ox/Ow)* + (Ov/Ow)? + (02/0w)? 
j0 on the lifting surface 99) Wwe have the wave equation, Eq. (19), in general coor T 
| at the leading edge dinates 
This singular part does not disturb Eqs. (26) and (27), _ Pa) ( U oO ) O ( V =) 

’ > T@ av > 70CeE 4 ’ : ‘se > “4 y ‘ , r 7 > , 5 i . » 
where we may replace K by AG. Phe singularity of ) ou \ VW du ov \UW ov 
the derivative G.’, giving the infinite part of the pres ; | 
sure p*, must be of the kind s *, if s denotes the dis- ° ( : ° )| mn Xa F=0 (33) su 
tance from the leading edge. The eigenvalue A may ow \U'V ow Ei 

a at: . ‘ > rq ve) , y yal i < oT rena° . . = ‘ : 
be calculated by inserting Eq. (28) into Eq. (27). his equation may be solved by separation of variables g1 

Both solutions G and G’ of the wave equation must : ; ‘ 

; aint : = 3 F = L,(u)-.M,(v)-N,(w) (34) 
vanish at infinitely great distances. Therefore, they 
must transform in the limiting case into the point solu- L,, J/,, and .V, are the wave functions. From Eqs. di 
tion (33) and (34) we obtain the differential equation 
r 
ten (Ut) M,,"(v) 4 WP V,,”(w) 4 UVW b(t) @ U P 
L,,(u) AM, (v) N,,(w) L,,(u) Ou \VIV va 
, , . % . res 
M,,'(v) =f J ) N,,'(w) O ( i ) 
— , —)}+x2=0 (35) an 
M,(v) Ov\ UN N,(w) Ow \UT 
All coefficients of this equation are known functions of u, if v and w have certain constant values, depending on the 
integral numbers ” and n’, if —n <n’ <n. Thus we obtain the function L,(u), and mutatis mutandis with the 
same constants the functions J/,(v) and J/,(w). 

The coordinate « may denote the family of finite closed surfaces beginning with the lifting surface 1 = 0 on the 
plane z = 0 and tending toward u = +. All orthogonal surfaces intersect the plane of symmetry z = 0 per 
pendicularly. Therefore, we have on this plane the derivatives 

: Ac 
0 Ou O . ra ‘ 
= . =. on the lifting surface 
Oz Oz Ou Ou 
(36) | 
Qo _w do _,, 2 — | rh 
= ° SS ° outside the hiting surface . 
Oz Os Ov ov ” lac 
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ny v=0 u=0 = Section Pressure distribution 
oth . 
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Fic. 2. Lifting surface with orthogonal coordinates 
51) 
The pair of wave functions 
S,(v, w) = M,(v) - N,(w) (37) 
32 denotes a surface harmonics of the finite closed surface « = constant. The surface element of the surface u = 
si constant is 
do = dvdw/VW (38) 
OT- waite ° ° 
rhus we have the orthogonality relation 
ca : 0, ifn 2 mi 
S,(v, Ww) + Sm(v, w)l(v, w)dv dw = J a (39 
E 11, ifm = m\ 
(v, w) denotes a function, which depends on the coordinate system. The integration extends over the whole 
33) surface « = constant. The wave functions 7, and N, must be normalized according to Eqs. (37) and (39). 
i For every integer m there are 2” + 1 different functions 1/7, and V,. We must choose odd functions 1V/,, V,, 
les giving opposite values at both sides of the lifting surface. A second restriction in the choice of the functions .1/,, 
V,, is given by the required kind of the singularity of G,, as mentioned above. 
mid rhe first kind of functions L,,(“) tends to infinity foru = ©. But a second independent kind of solution of the 
18. differential equation, Eq. (35), exists. We may write the Eq. (35) as 
= (35a 


L,"(u) + a(u) L,'(u) + b(u)L,(u) = 0 


Then we have also a solution!! 


> 


6 du’ ? ” ” 
R,(u) = L,(u) —— 9 4 — a(u" )du (40) 
Sus ie ) 


J u 


vanishing for“ = ©. We can determine the lower limit such that, for 1 ~ o, the limiting value, Eq. (30), 
This behavior is well known from the Bessel and Mathieu functions, the wave functions of the circular 


results. 
5) and elliptic cylinder. 
he (6) THe GENERAL FORMULA FOR THE KERNEL A 
ne , ‘ 
Now we are able to express the regular part of the kernel A, 

. R,,(u) 

1€ : j . = n ; ee 
G(u,v,v',w,w’) = >> ——— SAG, WISA0, (41) 

r- = 4 fee CO) 


According to Eqs. (36) and (41), the normal derivative at the lifting surface 1 = 0 becomes 
OG ; jO,ifv 2v’ orw2w’ | 
= U(0, v, w) > S,(v, w) S,(v', w’) = ie : r( (42 
(o,ify =v andw = w 


= 1 
This equation can be proved by means of the orthogonality relation, Eq. (39). The boundary of the lifting sur- 
face may be given by u = v = 0, see Fig. 2. At the front part (dotted) of this line, the function G,’ may have a 
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singularity, but not at the rear part. According to Eq. (42), we have the derivative for u = 0 
=i O,ifv20orwZ2w’ | 
G,, = U(O0~; 0, w) > S,(v, w)Syy(0, w’) = oe 7 os) 
: lo, ifv = Oandw u 


If we add symmetrical (s) and antisymmetrical (a) surface harmonics S,, we are able to cancel the infinity of G 
at trailing edge, required by the Kutta condition, see Fig. 2. This has been done by Schade’ for the circular lift 
ing surface. 

At the lifting surface we have u = 0, v < 0. 
6 function G, will not be disturbed, [compare Eqs. (42) and (43) ]. 


If we multiply G.,’ by a function /(v’) and add to G., then the 
Thus we obtain the complete kernel 


K = A[G(u, v, v’, w, w’) + h(v’)G,(u, v, 0, w, w’) | (44 
The plane z = 0 outside the lifting surface is given by v = 0, see Fig. 2. According to Eqs. (36) and (44), we have 
in this plane the derivative 
K, = AV(u, 0, w) [G,(u, 0, v,' w, w’) + h(v’)G,,-(u, 0, 0, w, w’) | (45) 
Further, we have 
x = x(u",0,w"); y = y(u’, 0, w”) = y(0, 0, w) 
ae On . : (46 
dx = du" + dw” = f(u", w) du 
Pa) ” alt 
cee 


where ~” and w” denote the integration variables. From Eqs. (24), (45), and (46), we obtain the function 


0 —(3 w/ 2) 
/ exp wx(u", 0, w)-f(u", w) V(u", 0, w)G,(u", 0, v’, w, w’)du” 
“ axe (47) 


7.7m 

J exp wx(u”, 0, w)-f(u", w) V(u", 0, w)Gyy(u", 0, 0, w, w’)du” 
- o + (¢ #/2 

The integration path 


(37) and 


where xv is the negative coordinate in front of the leading edge for real positive values of w. 
has been shifted in the complex number plane in order to obtain convergent integrals. According to Eqs. 
(41), we may write Eq. (47) as 


J F(u",w) > f.(u", v')N,(w)N,(w’)du” 
r 1 
(48) 


JS F(u", w) = 2(u")N»,(w)N»(w’)du” 


m ] 


hiv’) = — 


If we carry out the integrations and differentiate Eq. (48) with respect to w or w’, we obtain the infinite series 


= ) >, Hasm(v’, w, w’)*| 7 (49) 
Ow = gar ee N,(w), Nin(w) 


According to the Wronskian theory"! of the linear differential equations of the second order, we transform the de- 


N,,'(w), N,,'(w) Nn’ (w), Nm’ (wo) ” ” af (50) 
; . =I eee ‘exp | — a(w)dw 5 
N,(w), Nm(w) Nn(2o)s Nmm(20) ‘ 
the function a(w), see Eq. (35a). The determinant on the right-hand side of Eq. (50) vanishes at the axis of sym- 
wo, Where for antisymmetrical forms all functions V,, N,, vanish; for symmetrical forms, all their first 


terminant 


metry w = 
derivatives vanish. Therefore we have 


Oh/Ow = Oh/dOw’ = 0 (51) 


In order to obtain a simple formula, we put 


a2) = @) 


= Wo 
where wp») denotes the axis of symmetry y = 3 = 0 of the lifting surface. Thus we have 


or j Ox F du is 
x = x(u, 0, wo): dx. = du = — 52 
Ou U(u, 0, wo) ; 


wl 


Wi 


wh 


We 
for 


we | 
tion 


G 
lift 


the 
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47) 
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48) 
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From Eqs. (47) and (52), we obtain the final formula 


eo : . | (u, UO, w , 
exp wx(u, 0, Wo) G,(u, O, v', Wo, Wo) du 
; , - 2/2 U(u, 0, w 
h(v = “" : : 3 
iki ; V(u, 0, w 
exp wx(u, O, Wo)G,,(U, O, O, Wo, w du 
ees U'(u, 0, w 


This method is not restricted to the wave equation. It holds also for other linear differential equations of the 


elliptic type and second order in the subsonic range. But it may be more difficult to obtain the appropriate char 


acteristic function G for these differential equations. 
The kernel A has the difference property according to Eq. (12). Therefore, A should contain two independent 
sets of coordinates u, v7, w and uw’, v’, w’. But the characteristic function G(u, v, 7’, w, w’) does not contain wu’ 


Pherefore we have to modify the kernel A a little, while passing from Eq 


9) to Eq. (9a Following Eqs. (44) and 
17), we insert the product 
G,-(u, v, 0, w, w’)-G,(u", 0, v’, w, w’ 
into the integrand, uw” being the integration variable. Thus we must interchange uw and wu”; otherwise the kernel 
K’ Eq. (9a) would not satisfy the wave equation with the coordinates u, v’, w 
contains the variables u and u’ = 0. 


, 


It is shown below that the kernel of the pressure p*, given by Eq. (S83), 
Passing from Eq. (11) to Eq. (1la), we must also interchange u and wu’ in order to obtain a kernel that satisfies 


the wave equation. 


(7) THe Exvviptic LIFTING SURFACE” 


The orthogonal coordinates for a three-axial ellipsoid are the three roots p = X, uw, v of the cubic equation in p 
ix? (a? + p)] + [y?/(0? + 9) + [e*/(e + »)] — 1 = 0 5) 

where a* > 6? > c* are positive and 
—a@<cve —BPcpc—e<rA< 


\ = Xo yields an ellipsoid, 4 = wo, and v = v9 yield one-sheeted and two-sheeted hyperboloids. It is convenient to 


pute = 0. Then Eq. (54), with X = 0, vields both sides of the elliptic lifting surface with the boundary 


(x? a?) + (y?/b?) — 1 = 0; z= 0 (96) 
Further, we obtain 
. tf(r) : tf(u ; 1f(p) = 
Ul? = . : |? = s WV? = oe) 
(A — pw) (A — vp) (u — A) (uw — v (y—A)(v—yp 


f(A) = (a? + A) (0b? + A) (C2? + A) 


We may introduce other coordinates by the differential equations 


du = dd/2V f(A); dv = dp/2V f(u); dw = dv/2V fv »S 
which vield the elliptic function 9? of Weierstrass. Thus we have 
A= O(u) — s*; w= Ov) — s*; »v = O(w) — s*; 3s? = a? + O? 4+ Cc? 59) 


We introduce these coordinates into the wave equation, Eq. (35), and obtain the following differential equations 
for the wave functions: 
L,”(u) + [x??(u) + RO(u) + IIL,(u) = O j 
M,’(v) + [x??(v) + RYP(v) + 1)AM,(v) = 0 
N,”(w) + [k?Q?(w) + RP(w) + IIN,(w) = 0 \ 


60) 


k, / are numerical constants, depending on the integral numbers » and n’, if —n <n’ < n. By means of the dif- 


ferential equation of the Weierstrass function 
Y'(u) = 2V [P(u) + a? — s*] [P(u) + b? — s*] [P(u) +c? — 87] 61) 
we may express all (27 + 1) different wave functions L, as algebraic combinations of p(u) and a fundamental solu- 


tion of Eq. (60). 
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For the elliptic surface harmonics S,(v, w) the orthogonality relation, Eq. (39), becomes 


j0, ifn = m| 


i A S,(v, w)Sn(v, w) [P(v) — P(w) |dvdw = 1 ifn = 


(62 
= mf 


Instead of the function L,(“) we must use the second solution R,(u) according to Eq. (40). These functions 
must be inserted into Eq. (41) in order to obtain the characteristic function G and the kernel A. 
According to Eq. (54), there belong eight points (+x, +y, +2) to every set of orthogonal coordinates (A, u, v) or 


To avoid the inconvenience thereby produced, one may express the 7 function by the o functions of 
One may obtain G’ with- 


(u, v, W). 
Weierstrass.” With these functions the Kutta condition may be more easily satisfied. 
out the above-mentioned compensation procedure, as 1n the two-dimensional case. 


For incompressible fluids « = 0, the first term of the brackets of the Eqs. (60) vanishes. Thus we obtain the 
Lamé differential equation 
Ln’(u) + [ROP(u) + /][L,(u)] = 0 (63) 


where k = —n(n + 1). These potential functions L, were introduced by Lamé in 1837, while the wave function, 


Eq. (60), was studied first by Méglich'* in 1927 and M. T. O. Strutt® in 1932. 
(S) THE CrRCULAR LIFTING SURFACE 


For a = b, the elliptic functions reduce to the circular functions, 


O(u) = 1/cos?u; @(w) = 1/cos?v; O(w) = 0 (64) 
For the circular lifting surface, we have the following orthogonal coordinates: 
cos w sin w 
,=- _ y= - g =tgu-tghyv 
cos u-cosh v cos “-cosh v 
(65) 
cos” “-cosh v ; cos u-cosh? v , 
Uy ; Ve= - W = cos u-cosh v 
V sin? u + sinh? v V sin? uw + sinh? v 


We introduce these coordinates into the wave equation, Eq. (35), and obtain the following differential equations 


for the wave functions: 


Ln”(u) + [(x?/cost u) + (k/cos? u) + /] L,(u) = 0 (66) 
M,,"(v) — [(x?/cosh* v) + (k/cosh? v) + J]A,(v) = 0 (67) 
N,,"(w) + 1N,(w) = 0 (68) 


The differential equations, Eqs. (66) and (67), also define Lamé functions treated by Meixner.'** Eq. (68) yields 


N,(w) = exp iw V/1 (69) 


In this way the surface harmonics S,,(v, w) are simplified. 
For incompressible fluids x = 0, the first term in the brackets of Eqs. (66) and (67) vanishes. Thus the func- 


tions L, and A, reduce to spherical harmonics of the first and second kind. Schade'’ has treated this case in an- 


other way, but only for the six downwash distributions up to the second degree in x, y. 
(9) The INFINITE LONG RIBBON 
Another limiting case of the elliptic plate is the infinite long ribbon with varying downwash distribution in both 
coordinates x and y. We use the orthogonal coordinates 
x = —coshu-cosv; y= w; 2 = sinh w:sin vt (70 
r_9 ro : _— , ¢U) 
U-? = V-* = sinh? uw + sin’ 2; W=1 f 
Now we introduce these coordinates into the wave equation, Eq. (35), and obtain the following differential equa- 


tions for the wave functions: 


L,"(u) + [(k? — a?) cosh? u + /]] L,(u) = 0 (71) 
M,"(v) — [(x? — a?) cos? v + /] M,(v) = 0 (72) 


N’(w) + a?N(w) = 0 (73) 
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Eqs. (71) and (72) define Mathieu functions.'~* The constant / depends on n. Eq. (73) yields 
N(w) = exp taw (74) 
At first, we solve our problem for a special value of a and a corresponding downwash distribution. Let 
g(a) = (1/4) (x- — a (40) 
the usual parameter of the Mathieu functions. We obtain from Eq. (41) the regular solution 


. : : . ? > Ne,“ |q(a@), u | : - 
G(a: u,v, Vv, Ww, u exp 1alw a a. r <p S€n|G\ a), v}-se q\a),v ras) 
1 Ne,‘ ’[q(a@), 0] 


e, and .Ve, are Mathieu functions of the first and third kind.!'. The normalization is such that 


se,(0, v) = sin n 7 Vi 


It follows from Eqs. (70), for v = 0, 


¥/U = 1: X cosh u (7S 


/ exp(—w cosh u)G,(a; u, 0, 0’, w, w’) du 


/ exp (—w cosh “)G,,(a; u, 0, 0, w, w’) du 


: De ac . ; — 
A(a;u,v,v,W,Ww) = [G(a; u,v, v’, w, w’) + h(v')Gy(a; u, v, 0, w, w’) | (SO) 


Che exponential functions of w and w’ at Eq. (79) cancel each other. From this, by harmonic synthesis we ob- 
tain the sought-for three-dimensional kernel 


, 


, ; ; / = “ , 
A(u,v,v',w,w) = A > A (ma; u, v, v', W, Ww S1) 


which holds for any downwash distribution. The Mathieu functions for negative values of g(a) can be expressed 
by those for positive values of g(@).' This solution, Eq. (81), may be calculated numerically within the present 
knowledge of Mathieu functions. 


; 


If the given downwash w*(0, v’, w’) is different from zero only in a certain range —b < w’ < +4, then we have 
a rectangular lifting surface, oscillating between a gap of the infinite long rigid ribbon. In order to obtain the pres 
sure zero on the remaining parts of this infinite ribbon, one must add a certain opposite downwash distribution on 
these parts, which may be arrived at by an iteration method. 

It is interesting to note that, for the special values a = x, g = 0, Eq. (SO) reduces to the two-dimensional case 


+ for incompressible fluids, given below. 


(10) THE Two-DIMENSIONAL CASE 


By putting a = 0 in Eqs. (75) to (80), we obtain the two-dimensional case of the lifting surface theory. Thus 
we have A = 2/7, and we may write the velocity potential according to Eqs. (9) and (75) to (SO), 
2 
o*(u, v) = — £5 w*(0, v’) sin v’ dv’{ G(u, v, v’) — 
t us 
/ exp(—w cosh u) G,(u, 0, u’) du 
. y + (4 #/2 
G,’(u, v, O) = ol (S2) 
1 « 
} / exp (—w cosh u) G,,’(u, 0, 0 du 
} e + ta/2 


| For the first integral, Cauchy's principal value has been taken. By applying the operator |w + (0 Ox)] to ¢* and 
using known properties of the Mathieu functions, we obtain the pressure according to Eq. (11 
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2 ‘ , ° / / j ) ¥ / 
p*(u,v) = — w*(0, v’) sin v’ dv oo — ——— —_ Ga, 0, 2") + 
wT. sin v’ Ov 
{ 


G,(u, v, 0) [G,(0, 7, v’) T(x, w) — G,(0, 0, v’) | (83) 


(t w/? 
/ exp (—w cosh u“) G,,’(u, m, 0) du 
Se TE (84) 


T(x, w) = = 7 a7 
/ exp(—w cosh u) G,,’(u, 0, 0) du 
F » + een 
The author" has proved that Eqs. (82) to (84) are identical with those of Timman and van de Vooren" and of E. 
Reissner.'® These authors obtained their solutions in another way and in a complicated form, being valid only for 
the pressure on the lifting surface u = 0. 
For incompressible fluids x = 0, the characteristic function G reduces to 
, F l cosh u — cos(v — v’) 
G(u,v,v) = n - 
4 cosh u — cos(v + v’) 
Thus the Eqs. (83) and (84) transform term by term into the known solution of Kuessner and Schwarz." 
The velocity potential ¢* according to Eq. (S82) has not been published until now, nor has the case « = 0, al- 
though it may be useful for the airplane stability theory. Inserting Eq. (85) into Eq. (82), we obtain the ve- 


locity potential for x = 0 


‘ I —_" ee ee |, cosh “ — cos (v — v’) 
o*(u, v) = w*(0, v’) sin v’ dv lt rl 
0 


1 
va (2 cosh u — cos (v + 2’) 
© — (4/2) : du 
: Fes, exp(—w cosh u) , 
sin v sin v © + (ix/2) cos v’ — cosh u sa 
ee - (86) 
cosh “u — cos v e du 
exp (—w cosh 1) 
© + (ix/2) 1 — cosh u 
According to Watson, '* we have the formulas ) 
sa (iw/2) 
‘ —y) 1 > ° 
/ exp(—w cosh u) cosh mu du = mi” HT, ( ~ i) | 
2+ (tr/2 
(87) 
0 — (in /2 
exp (—w cosh u) sinh mu du = 0 
0» + (tr /2) 
From this we obtain the formula 
eS eree , du ; ae er ae 
exp (—w cosh 1) = twr [+ Ho? (—tw) + 1 (—tw) | (88 
» + (ix/2) 1 + cosh u i 
which we must insert into the Eqs. (84) and (86) for « = 0. 
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 e. 
m= ‘ 4 Subscripts 
for ABSTRACT 
h = 1n translational motion 
Semiempirical results by Halfman, et al.,? are employed to J = imaginary part 
demonstrate the existence of steady-state oscillations (limit j = 90 deg. out-of-phase component 
cycles) of pitching airfoils in separated flow (stall flutter). The r = characteristic of dynamic state 
gn flow around an oscillating two dimensional airfoil was studied R = real part 
with a Mach-Zehnder interferometer. This investigation shows a = in pitching motion 
the prime effects of frequency k = (wh)/v, initial angle of attack 
a;, cascade effect, and Mach Number on the aerodynamic reac 
al tions in separated flow INTRODUCTION 
al- | 
ve- Ballenger VIBRATIONS OF AIRFUILS at high 
SYMBOLS . i ‘ Mag 
angles of attack were first explored by Studer’ in 
b = semichord of airfoil 1936. Studer recognized that the separation effects at 
CL = L/4gb, lift coefficient stall were involved, but the impact of the successful clas 
E = WM /4gb?, moment coefficient ‘ 
” : eae sical flutter theory was -uch as to prolong the attempts 
‘ = 2b, chord of airfoil ; s : 
h = translational deflection of center of torsion, positive to correlate the experimental data to . slightly change¢ 
downward classical theory. <A detailed study of the real cause of 
i = V-1] the phenomena was published by Halfman, et al.,* in 
86) k = (wh)/v, dimensionless frequency 1951. They investigated the aerodynamic derivatives 
I - aerodynamic lift in center of torsion per unit span at high angles of attack. It was clearly shown by 
of airfoil, positive downward ; . . eae 
; strain-gage recording of lift and moment how the sep 
| V/ = aerodynamic moment about center of torsion per : : : 
j unit span of airfoil, positive when working in the aration changed the aerodynamic Feactiouls at stall, 
same direction as a and a semiempirical technique was devised for the pre 
VW = Mach Number diction of the aerodynamic moment loop of pitching 
) = y2/2) “te ic ressur . ae 2 : 
q = p(v?/2), dynamic pressure motions in the separated region. 
87) Re = Reynolds Number 
= distance between two blades in a cascade In the case of turbomachine blading, the author has 
T; = characteristic time shown that the classical flutter is unlikely to occur, 
. = UUme since the aerodynamic forces are too small to change 
= free-stream velocity : ‘ : : 
; ; eS the normal modes into a combined pitching-transla 
Ware = work per cycle due to moment, positive when going ; . , 3 
sto the tiade tional motion at a common frequency. Compressor 
88) Wrre = work per cycle due to lift, positive when going into and turbine blades, therefore, will vibrate in their nor 
the blade mal modes with their natural frequencies very little 
Xp = chordwise distance from leading edge of center of changed. The vibrations will eventually reach a 
torsion of airfoil _ 
| ; a steady state (though the related stresses of course may 
a = angle of attack, counted from angle of static equilib : ah F 
rium a, exceed the fatigue strength). 
a; = initial angle of attack at static equilibrium (flow This indicates the existence of limit cycles, in clear 
direction-chord direction ) _ ‘ . : We - ‘ 
‘ distinction from the arbitrary conditions of state in the 
ers a = amplitude of angle of attack (a; + apo) ‘ : 
6, A = slope of quasi-steady moment coefficient classical flutter where only the borderline case between 
é = stagger angle in cascade (chord direction axis through a damped and an unconditionally divergent state may 
o- 7 similar points in cascade ) be obtained from the linear system. The existence of 
pp. p = eres — ies these limit cycles was first discussed by the author in a 
o = phase shift of moment coefficient ai 4 ¢ ‘ ‘ = 4 , " 7 ; 
7 ; preliminary report in April, 1952, at the Gas Turbine 
w = frequency, rad. per sec . iti atte . : ; 
ch, Laboratory at M.I.T. The experimental investiga- 
— Received June 1, 1953 tion reported in the present paper demonstrates the in- 
* . . a  . p e . “ted : Cac ” e : e ° 
} _* Abstract of part of an investigation conducted at the Gas Tur- — fyence of various factors on the separation phenomena 
tial | bine Laboratory of M.I.T. in partial fulfillment of the require . ; . ain Nn a 
5. and the forming of lift and moment traces which indi- 
can ments for the Sc.D. degree : ; : ‘ 
cates that the flow feeds mechanical work into the vi- 


t Assistant Professor of Aeronautical Engineering. Presently 
with the STAL Company, Finspong, Sweden. brating blade system in the separated region. 
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Fic. 2. Calculated hysteresis loop of aerodynamic moment 
in stalled flow. Initial angle of attack, a; = 18°; amplitude, 
a = 10°; and frequency, ka = 0.10. 
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Fic. 3. Calculated hysteresis loop of aerodynamic moment 
in stalled flow. Initial angle of attack, a; = 18°; amplitude, 
a = 10°; and frequency, ka = 1.0. 
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(1) THe CONCEPT OF THE Limit CYCLE 


Consider a two-dimensional airfoil elastically re 
strained so that either a pitching motion around a fixed 
center xp or a pure translational motion may take place. 

The work per cycle transferred from the flow to the 
airfoil is determined by the following integrals: 


Wur = f M dea (1 
Wir = f Ldh (2 


These expressions represent work fed inte the vibrat- 
ing blade when a closed area in the .W/-a or the L-h- 
plane is circumscribed in a clockwise sense in the course 
of a cycle. 

In order to establish the form of the moment trace 
in the A/-a-plane (pitching motions), the Halfman 
group’ investigated the classical potential flow case 
and then transferred or modified the results to the case 
of separated flow. Their hypothesis was that the value 
of Cy at an instantaneous angle of attack differs from 
the static-state value at the same angle of attack by an 
amount ACy. It was assumed that the rate of return 
of Cy toward the static state is proportional to the 
instantaneous difference AC, (exponential approach). 
Therefore, 


~ 


dCy dt = —~ACy i, (3) 
where 7’, is a time constant of the process. 

In this hypothesis, however, the static-state values 
were not the usually measured values of the moment 
curve Cy = Cy(a@) but were represented by a dynamic 

The difference ACy was referred 
This line intersects the true static 


line with the slope X. 
to this line, Fig. 1. 
curve at the mean or so-called initial angle of attack 


1 — i(k/k;) ~ 
- d+ aye 1) 
1 + (k/k,)? 


Q; 
It was found that, for potential flow, 


Cu — Cu; = real | 


The parameters \ and k, are functions of the fre- 
quency k = wb/v of oscillation. Eq. (4) corresponds to 
an elliptic trace in the 1/-a-plane (see Fig. 1). This 
elliptic trace encloses an area that, according to the 
conditions of Eq. (1), represents a negative work—that 
is, the aerodynamic moment exerts a damping on the 
pitching blade. 

When the instantaneous angle of attack is larger than 
the angle of stall at static condition, it was found in 
reference 2 that the moment ceases to follow the elliptic 
path, a more or less precipitous change of the moment 
takes place, and an area in the .\/-a-plane is enclosed 
in the opposite sense, indicating that the aerodynamic 
moment now has fed energy fo the vibrating blade (see 
Figs. 2 and 3). 

The same general hypothesis of an exponential ap- 
proach to a ‘‘steady”’ state was used in the stall region. 
However, the line of slope \ was replaced by a number 
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FLUTTER OF 


‘ straight lines chosen in a manner that made a com- 
yuted cycle to follow the experimental findings as 
closely as possible. The position in the cycle where 
the moment trace starts to deviate from the elliptic 
path of the potential flow as a result of separation must 
be obtained from experiment, because there is a certain 
lag that makes the ‘‘dynamic’”’ stalling angle larger 
than that of the static condition. 

The above method was used by the author despite 
the tentative nature of the assumptions, the main thing 
being that the same ‘‘guiding’’ lines were used consist- 
ently in the stall region to replace the A-line. The 
“blunt” airfoil of reference 2 was the object of investi- 
gation in the following study of the limit cycles. It 
had a static moment coefficient Cy rising to 0.126 at 
a; = 14° and falling to zero between 15° to 18°. 

Some previous work* ° showed that the observed stall 
flutter of compressor blades appeared to reach some 
steady-state harmonic vibration. This suggested that 
it would be possible to calculate the amplitude of these 
vibrations. The amplitude of such limit cycles is, of 
course, defined by the criterion 


Wur = { Mda=0 (6 
Wir = f Ldh=0O (6) 


The physical existence of these limit cycles may be 
explained thus. At stall, the blade picks up a disturb- 
ance that feeds work into it and so increases the ampli- 
tude of vibration until a state is reached where the blade 
sweeps over a sufficient range of nonseparated flow ex- 
erting damping on the blade. The balance of Eqs. (5) 
and (6) is then fulfilled. 

Since the initial angle of attack a; and the natural 
frequency k, appear to be of prime importance, this 
suggested that the simplest relationship of the steady- 
sta'e amplitude a) (for constant Mach Number and 
Reynolds Number and unchanged blade geometry) 
would be written 


ay = ao(aj, Ra) (7) 


The general result of the Halfman group indicates 
that the effect of stall will culminate for a certain region 
of a,- and k,-values. This may again be understood 
from a simple physical consideration. 

When a; is very different from the stalling, then the 
airfoil is either not stalled at all or completely stalled in 
its initial state. In both cases, it takes a considerable 
angular deflection (connected with the usual damping 
reactions) to bring about that change of the aerody- 
namic moment at stall which is responsible for the 
exciting action. Because such larger deflections are 
connected with considerable aerodynamic damping, it 
is a fair conclusion that beyond a certain a;,-range the 
exciting action at stall is not able to sustain a deviation 
from the static state large enough to retain the con- 
tact. 
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Fic. 4. Calculated variation of steady-state limit cycle ampli 
tude with variation of frequency k. Pure pitching motion, x, = 
37 per cent chord and Re = 1 X 10° 


Concerning the frequency k,, there exists an upper 
limit (which may vary with a;) above which the bound- 
ary-layer separation is slow in comparison with the 
motion of the airfoil and no separation has time to 
develop. On the other hand, for small k,, all devi- 
ations from the semisteady state are small and no 
large loops are formed at all. The enclosed area in the 
M-a-plane, representing work per cycle, then becomes 


small. 

A number of loops were calculated for a; = 18 
and 14°. In order to find the limit cycles defined by 
the steady-state amplitude aj at which Wyre = 0, 


it was necessary to calculate complete loops with both 
negative and positive work per cycle for any given pair 
of k,- and a;-values. 

Examples of the obtained moment loops are shown 
in Figs. 2 and 3 for the values a; = 18°, a = 10°, and 
k, = 0.10 and for a; = 18°, ag = 10°, and k, = 1.0. 

The important relation between the steady-state 
amplitude ag and the natural dimensionless frequency 
k, = (wgb)/v is reproduced in Fig. 4. The spread or 
uncertainty of the calculations is indicated by the 
shaded region. The curve is shown for a, = 18°, and 
one point for a; = 14° is also given. 

For each angle of attack a; in the vicinity of stall 
there exists a similar curve. Taken all together these 
curves form the relation ap = ao(a;, Ra). 

It requires only a few simple operations to perform a 
harmonic analysis in order to calculate the equivalent 
harmonic components of Cy in phase with a and da/dt, 
when the usual moment loop in the 1/-a-plane is known. 
Denoting the component of .1/ in phase with the deflec- 
tion ./, and that in phase with the velocity da/dt, M,, 


it is found? that 
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The commutator system was later 


Fic. 5. 


View of test section. 
mounted on the shaft to the right in the picture 
M, = Wua/e-a0 = $M da/x- a (8) 
Mr = Was/r- a = far da;/m- ao (9) 
Here, 
a = ao'cos wl (10) 
aj; = ao'sin wl (11) 


A few computations of the Cy-components were 
made at a; = 18° and k, = 0.3 (the peak condition in 
Fig. 4). Since it was found from reference 2 that the 
separation (incidence of dynamic stall) could take place 
within a certain region of ‘“‘overshoot,’’ two different 
values of this overshoot were investigated at ag = 16 


The results are given in Table 1. 


TABLE | 
Calculated Cu-Components in the Vicinity of the Steady-State 
Condition at aj = 18° and ka = 0.3 
a CMR CMI V Cmr? + Cui? OM War 

16° (early stall) 0.0525 0.0648 0.084 51 
16° (stall at maxi- 0.155 0.013 0.155 5 

mum angle of at- 

tack) 
24 0.1667 0.011 0.167 -3.8 - 


The method of representation in Table | coincides 
somewhat with the way in which the condition of flutter 
has been discussed in the past. It shows how the prob- 
lem of dynamic stability no longer is to determine the 
discrete k,- and v-values at which an unlimited vibra- 
tion may start. Instead it is found that, for any given 
pair of a,- and k,-values, there exists at least one value 
of the amplitude ap (between 16° and 24° in the above 
example) for which the phase shift ¢, = 0. At this 
condition, the equivalent Cy-vector is exactly in phase 
with the motion vector a, and, therefore, no energy 
(J Ma dt = 0) is transferred to or from the sys- 


tem. 
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A general investigation of the dynamic stability of 
single degree of freedom systems based on a simplified 
behavior of the reactions at stall has been made by the 
simulated by stepwise 
These jumps occur 


author.* The separation is 
jumps in the lift and the moment. 
in general in different parts of the cycle, depending on 
whether the angle of attack is increasing or decreasing, 
thus allowing for the hysteresis effect. 

It is thus possible to show that many flutter motions 
require an initial impulse of kinetic or potential energy 
to an elastically restrained blade in order to start 
(hard flutter). Regardless of the initial state, however, 
Fig. 4 explains why reports on compressor blade flutter 
usually set a value of the measured stress (at the root 
section, for example) which defines a critical k-value.* 
This is equivalent to defining an arbitrary limit ampli- 
tude ap (horizontal line in Fig. +) which intersects the 
amplitude curve at the “‘critical’’ value. An upper 
limit of k, beyond which no appreciable change in the 
degree of separation is likely to occur at all, is discussed 
in the following section on the experiments. 


(Il) EXPERIMENTAL RESULTS 


Test Setup 

The test section was mounted in the cascade branch 
of the high-speed, variable-density wind tunnel at the 
Gas Turbine Laboratory of M.I.T. (Fig. 5). 

The cascade consisted of five blades of NACA 6.5409 
sections mounted with a spacing s/c = 1.0 and a stagger 
angle € = 52° between the chord direction and the 
direction of a line through similar points of the cas- 
cade. 

The angle of attack of the four extreme blades was 
fixed to 13°, whereas the middle blade could be set at 
any incidence on the supporting drill rods (Fig. 5). 
These drill rods could be oscillated mechanically in a 
steady harmonic motion by an electric motor system. 


Both translational and pitching motions could be 
obtained. A 5-in. Mach-Zehnder Interferometer was 


used for observation of the flow, and single exposures 
were obtained in connection with a synchronization 
system (commutator on the electric motors) which 
triggered a high-intensity spark-gap-type light source 
with magnesium electrodes. 

The interferometer gives a possibility to obtain lines 
of constant density and, from these, lines of constant 
pressure in isentropic flow. The density lines are either 
directly obtained by the method of infinite band spac- 
ing (Fig. 9 and upper exposure of Fig. 10) or by a super- 
position of a flow and a no-flow picture. This latter 
method is demonstrated in Fig. 11, which shows the 
flow picture and the superimposed picture. 

Eight exposures were taken in each run (closed loops 
in Fig. 7), and the lift and moment traces were faired 
in between the calculated values for these eight posi- 
tions. 
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RESULTS 


Some results are shown in Figs. 6 and 7 of the static 
and dynamic lift and moment characteristics in pitch- 
ing motion. Fig. 8 indicates the growth of the bound- 
ary layer or wake in one kind of pitching motion. 
Figs. 9 and 10 are interferograms of pitching motions, 
and Fig. 11 shows an exposure of a translational mo- 
tion. Some of the main findings will be discussed be- 


low. 


Effect of Frequency k = wh/v (Pitching Motion) 


This effect was the subject of discussion in Section 
(1). It was found that small k-values should corre- 
spond to cases with the dynamic characteristics in 
close agreement with the static lift and moment curve, 
both in the unstalled and stalled regions. With in- 
creasing k, there would be an increasing deviation from 
the static behavior. For still larger k, it is to be ex- 
pected that this deviation corresponds to a suppression 
of the separation effects in the stall region (the time 
constant of the boundary-layer separation large in com- 
parison to the time 7 of one cycle of oscillation) and 
a tendency toward the classical elliptical trace of Figs. 
1, 2, and 3 over the full course of a cycle of oscillation. 
The actual frequency k of the experiments was small 
(to the left of the resonance peak of Fig. 4). Accord- 
ingly, we may expect that the separation effects will 
increase with increasing k, because we are in a range 
of k where the separation effects are delayed but not yet 
suppressed. 

Before a comparison can be made, it is proper to 
point out that in the runs in cascade (s/c = 1.0) there 
was a lower-surface separation of the same kind as 
for the upper surface (compare the upper part of Fig. 
9 for the lower-surface and the upper part of Fig. 10 


for upper-surface separation). 


One run was made at k = 0.023, 17 = 0.57; another 
was made at k = 0.029, \J = 0.48; and a third was 
made at k = 0.08, \J = 0.27. The mean angle of 
attack a; for these runs varied between a; = 11.5 


to 14.5°, and the center of torsion x, was at around 42 


per cent of the chord. 

Studying the dynamic characteristics of Fig. 7, one 
will find how the delaying effect of the motion brings 
about increasing hysteresis loops as k increases from 
0.023 over 0.029 to 0.08. This is true for both lift and 
moment. It means that for these k-values the real 
steady-state limit amplitude a must be larger than the 
experimental amplitude. 

A study of the boundary layer for the case a; = 
13.7°, ao = 12.9°, and k = 0.08 with WV = 0.27 is of 
special interest because it indicates the effect of motion 
(Fig. 8). 

The phase lag of the separation in the dynamic case 
relative to the steady state was found to be of the order 
(1/9)-m rad. (27 rad. per cycle). Assuming that the 
separation occurs at the maximum angle of attack, a 
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Fic. 6. Static lift and moment characteristics of pitching mo 


tions. (NoTE: The neighbor blades were fixed at a; = 13° at 
these measurements 


phase shift of more than 7 rad. corresponds to a case 
where the reattachment of the boundary layer takes 
place at increasing angles of attack. It is therefore rea 
sonable toexpect that a phase lag of 7 rad. (one-half cycle 
of vibration) represents a limit beyond which no appre 
ciable separation occurs. The corresponding limiting 


0.08 /(1.'9 0.72. It is inter 


frequency is (R,) 
esting to compare this value with Shannon's report‘ 
that the average frequency at which flutter starts is 
k = 0.75 when starting is indicated by measured stresses 
equal to +2 tons per sq.in. in the 1- by 0.5-in. cantilever 
blade under investigation. This means that an investi 
gation of the time necessary for boundary-layer separa 
tion on one hand and the practical case of a vibrating 
cantilever blade on the other both indicate the order of 
magnitude of the frequency k, beyond which the whole 
separation effect will not occur at all. It is likely that 
this limit is rather vague—that is, there is a region of 
k-values in which separation occurs sporadically, even- 
tually with random distribution in time. 


Influence of Initial Angle of Attack (Pitching Motion) 


This factor was investigated by a comparison of two 
runs at k = 0.029 and J = 0.48 where the initial angle 
of attack of the middle blade was a; = 11.3° and a; = 
21°, respectively, whereas the fixed neighbor blades 
had a; = 13.0° in both runs. The first run at a, = 
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Fic. 7. Dynamic lift and moment of pitching motions in cascade. (Note: The neighbor blades were fixed at aj = 13° at 


these measurements. } 
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UPPER SURFACE Phis means that one may expect a beneficial effect 
by making the spacing closer. 











0 17/3 2/37 

+2 i There was some indication of the choking effect, item 
3° 7 , ; 
>» 40 ‘a 2 (3) above, and this was particularly noticeable in the 
a te -a! | translational runs. Thus, the lower surface stall was 
2 u | ° k | . ° ° op e : , 
x . e canceled in the case of a single airfoil with removed 
=< 20 os : ors . . 
x | 75%CHORD| & = 0.08 neighbor blades. There was also in some of the trans 
mi AT | lational runs some trend toward a somewhat larger 
us “i | 

| a 
3 oO 
o | 
x 
< 75 % CHORD| 
S 
5 
°o 
oO 
7 4/37 5/37 27 


LOWER SURFACE 
Fic. 8. Wake thickness of upper and lower surface at 75 per 


cent chord distance. Pitching motion kg = 0 and 0.08, a = 
13.5°, a9 = 13°,and M = 0.27. 


11.3 is identical with the middle run in the discussion of 
k, whereas the run with a; = 21.0 is distinguishable to 
the outermost right in Fig. 7. It is seen that despite 
the constant frequency there is a difference in the 
change in lift and moment in the stall region. This is 
naturally found to arise from a difference in quality 
of the stall. Usually the separation from the upper 
surface is moderate, and the interferograms show only 
a slight displacement (see upper part of Fig. 10, and 
compare the order of magnitude of the boundary-layer 
thickness 6/c in Fig. 8). 

For the larger angles of attack, however, it was found 
that the whole region above and beyond the upper sur- 
face was vigorously disturbed. The general flow pic- 
ture was like that of the lower part of Fig. 10, which 
shows a single airfoil for which, however, the same 
amount of separation occurs at an angle about 10 
smaller than that of the cascade (compare the static 


characteristics of Fig. 6). 


Cascade Effect 

The influence of the cascade on the actual blade may 
occur in at least four different ways—namely (1) effect 
of static state, nonuniform velocity field due to static- 
state circulation of the blades in a cascade; (2) effect 
of the vorticity along the vibrating neighbor blades 
and their wakes (each blade performing a small har- 
monic motion); (3) choking effect when the channel 
between two blades is periodically diminished; and 
(4) the transient effect of stall from blade to blade 
which causes traveling stall along the cascade (rotating 
stall propagation). It follows from the description of 
the test arrangement that not all of these effects could 





be investigated. 
= , . ‘ Fic. 9. Interferograms of pitching motion in cascade, s/c = 
First of all, it was found that the cascade in general 1.0, (a; = 13.0°, a = 13°, ka = 0.029, and MV = 0.48 Upper 
has a suppressing effect on the stall (compare static exposure: ay + a = | o, increasing incidence i vue re xp Sure 
. ai + a = 23.3°, decreasing incidence. Flow direction parallel to 


characteristics of single airfoil and cascade in Fig. 6). the vertical wire. Infinite band spacing at no flow condition 
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Fic. 11. Interferograms of translational motion in cascade s/< ca 
= 1.0. Upper exposure: a; = 26.3° and M = 0.62. The blade 
is displaced downward 6 per cent and moving toward its neutral 
position with the frequency k,= wbh/v = 0.016. Flow direction m1 
parallel to the vertical wire. (Note the shock wave at the 
leading edge.) Lower exposure: The same picture as above but 
superimposed on the no-flow fringe pattern to obtain the band Wi 
shift of constant density (blurred lines). This technique was 
used in principle but improved in order to calculate the density 
and related pressure distribution over the airfoil. se 


la 
b 





Fic. 10. Interferograms of pitching motion. Upper exposure: 
(a; = 18°, ao = 18°, ka = 0.029, and M = 0.48, cascade s/c = 
1.0), a¢ + @ = 28.2°, maximum incidence. Infinite band spac- 
ing at no flow condition. Lower exposure: (aj = 15.8°, ao= 13°, 
ka = 0, and M = 0.48, single airfoil), aj + a = 28.8°, maxi- 
mum incidence. Finite band spacing at no flow condition. Note 
the wide, disturbed wake over the upper surface 


wa. 
— 


% er 


¥t 
5) 


A 





5 td 
lade 
tral 
tion 
the 
but 
and 
was 
sity 


FLUTTER OF 


separation effect when two blades were closer together. 
his was more pronounced at the higher Mach Num- 
bers. 

The first type of cascade effect, item (1) above, was 
suggested by Bellenot and d’Epinay’® to be a first-order 
cause, which alone should be enough to explain the 
single degree of bending flutter. In the present experi- 
ments, however, it was found that, as long as the sep- 
aration was of the same type (that is, equally large 
wakes), there was a rather small influence on the lift 
caused by the varying position in the cascade. This 
seems to confirm extended calculations by the author, 
which showed that no strong influence could be ex- 
pected. Therefore, the experiments of the author, as 
well as recent experiments on the rotating stall propa- 
gation, seem to indicate that the importance of the cas- 
cade is primarily connected with the effect the cascade 
has on the separation phenomena on the surface of the 
airfoil. 

For a single airfoil, it is known that only the velocity 
of the translational deflection—that is / (or rather the 
but not the de- 


apparent angle of attack a’ = h/v) 
flection itself controls the aerodynamic forces. If there 
is an influence from the part of the cascade on the stall 
(choking effect, for example), then it is apparent that 
the position also has an influence on the aerodynamic 
forces. 

Now a separation effect due to the influence of a’ = h/v 
alone is enough to explain the possibility of single de 
gree of freedom flutter. The present experiments defi- 
nitely confirmed that, when the airfoil is working all 
the time in translational motion in separated flow, the 
aerodynamic lift will be able to feed energy into the 
The influence of the cascade may be some sort 
In order to find this influ- 


blade. 
of control over the stall. 
ence, one has to make more extended investigations. 
The order of magnitude of the forces involved in trans- 
lational motion is roughly indicated by the difference 
between the maximum and minimum measured values 
of lift in the full course of a cycle of vibration. This 
difference Cymer. — Crmin, WaS found to be of the order 
of 0.10 in the lift coefficient of translational motions 
with h/c < 0.23 and represents essentially different 


degrees of separation during the harmonic cycle. 


Effect of Mach Number 


This effeet has often been rather overlooked in the 
past research of turboblade vibrations. The general 
effect appears to be a decrease of the maximum value 
of the lift coefficient C, and moment coefficient Cy 
(compare static values in Fig. 6 for pitching mo- 
tions). 

The difference Cryo, — Crmin. aS discussed above had 
a slight trend to culminate for medium Mach Numbers, 
the small Mach Numbers showing a more stable flow 
with smaller effect of stall and the larger subsonic Mach 
Number (0.6—0.7) corresponding to a separated, but 
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more consistently separated, flow with smaller internal 
changes (Cymer. — Cimin. ™ 0.08). 

The shock-wave pattern visible in Fig. 11 at the lead 
ing edge points toward important Mach Number ef 
fects (compressibility burble connected with so-called 
shock flutter). The type of shock usually attributed 
to the laminar boundary layer (A-shock), as well as the 
short strong shock of turbulent boundary layer, was 
found. However, that does not necessarily mean the 
same thing as usual in this particular case of transient 
boundary layers. 

The importance of these phenomena justifies a sep 
arate investigation of somewhat higher Mach Number 
and smaller angles of attack with a continuous recording 
of the flow pattern, with a motion-picture camera, in 
order to find an accurate relationship among the shock 
pattern, the boundary-layer separation, and the mo 
tion of the airfoil. Now it was found that the boundary 
layer always was separated behind the shock waves, 
but it was found to be difficult to correlate any differ 
ence in the amount of flow separation with the type 


of shock wave observed. 


Effect of Shape of Blade Profile 


Though only a single profile shape was tested, there 
was an indication of the possible influence of the blade 
profile, given by the appearance of lower- and upper 
surface separation. The sharper profile form of the 
lower surface around the leading edge caused sepa 
ration in cascade at an angle of attack around zero deg. 
(velocity direction parallel to chord). The turbulence 
was of a medium level (sphere drag coefficient had 
Re-ric, = 195,000). The Reynolds Number of the blade 
was Re = w-c/v = 500,000 at MV = 0.5. 


CONCLUSIONS 


(1) The existence of steady-state limit cycles in pitch 
ing motion at stall has been shown by use of a semi 
empirical method for calculation of the moment loops. 
The result in its simplest form is expressed by the rela 
tion ap = aola;, k. ‘ 

(2) An interferometer study confirms the existence of 
exciting aerodynamic reactions in pitching and trans 
lational motion at stall. The extension of the sepa 
ration phenomena is shown to be influenced by such 
variables as k, a;, \J, cascade geometry, and airfoil 
shape (the effects of Re, turbulence, and vibrating neigh 
bor blades were not studied). 

(3) A study of the growth of the wake over the upper 
surface on a blade in motion indicates the order of mag 
nitude of the time constant of the separation and ex 
plains the existence of an upper limit of k, (and most 
likely of k,), above which the airfoil motion is too fast 
to allow any separation phenomena to take place, which 
seems to be in agreement with experiments of vibrating 
cantilever blades.‘ 


(Concluded on next page) 
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‘Temperature Distribution in Laminar 
Stagnation-Point Flow with Axisvmmetryv 


CHIA-SHUN YIH* 
Towa Institute of Hydraulic Research, State University of Towa 


ABSTRACI 


Based on the velocity distribution given by Homann in 1936 
for the laminar incompressible flow near a stagnation point in 
axisymmetric flow, solutions for the following cases of forced 


convection are obtained: 

(a) A point source of heat is situated at the stagnation point 
on an insulated boundary 

(b) The difference between the temperature of the boundary 
and the ambient temperature varies as an arbitrary power of the 


distance from the stagnation point 
(c) The temperature rise is a result of viscous dissipation 


alone 

Solutions for practical cases are obtained by combining the re 
sult for case (c) with that for case (a) or (b 

Numerical results are given for all three cases for various para 
metric values of the Prandtl] Number, with special emphasis on 
the Prandtl Number 0.7, which applies approximately to air 


under normal conditions 


(1) INTRODUCTION 


o bw VELOCITY DISTRIBUTION in laminar stagnation 
flow with axisymmetry was obtained by Homann! 
in 1936. Since his solution is one of the rare exact 
solutions of the Navier-Stokes equations of motion, it 
seems desirable to study the temperature distribution 
in the fluid for various boundary conditions—on the 
assumption that the velocity distribution given by 
Homann is essentially undisturbed by the nonhomo- 
geneous temperature of the fluid. 

With the velocity distribution determined from the 
Homann solution, the energy equation can be solved 
for various boundary conditions. Temperature rise 
due to viscous shear can be discussed separately. 


(2) VELocITY DISTRIBUTION 


Since the velocity distribution obtained by Homann 
is to be used in the calculations, a brief résumé of his 
solution is necessary. If the stagnation point on the 
bounding plane is taken as the origin and if 7 and 2 are, 


respectively, the radial distance from the origin and the 


vertical distance from the boundary (see Fig. 1), then 
the Navier-Stokes equations of motion are 


Ou Ou 1 Op 
u + w =_— 
or Oz p Or 
(= 1 Ou u ~~) ' 
) _ (1) 
or? r Or r? 2? 
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Ow ze) 1 Op 
u + w = << + 
or Oz p Oz 
O*w 1 Ow O-w 
y(— + + ) (2) 
or- r Or Oz? 


in which uw and w are the velocity components in the r- 
and z-directions, respectively; p and y are, respectively, 
the density and the kinematic viscosity of the fluid; and 
p is the pressure. The equation of continuity is 


(O Or) (ru) + (0 Oz) (rw) = O 


which permits the use of the Stokes stream function ¥ 


such that 


| oy 1 dy 


“ = WwW = (-) 


r Os’ r or 


If one takes 


n = (28/»r) “’s | 

in which @ is a constant, one obtains from Eq. (3 

“u= Brf"(n), w= —(28vr) ‘f(y () 
where the prime denotes differentiation with respect to 
n. Since 8 is so far not specified, one can take f’( 
1. The radial velocity just outside of the boundary 
layer is then “7, = 67, and 
l Op On; 


i; = pr a) 


p Or Or 


It then follows, from Eq. (1), that 
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Fic. 2. 


f'2 — Off" = 1 + af" (7) 


which is to be solved with proper boundary conditions. 
After f is obtained, one can integrate Eqs. (2) and (6) 
and obtain an expression for the pressure. The equa- 
tion obtained and solved by Homann differs from Eq. 
(7) only in that the coefficient of f’’’ is 1 in Homann’s 
equation instead of 2, because of the absence of the con- 
stant 2 in the expressions assumed by Homann for y 
and y. The foregoing slightly modified treatment was 
used by Fréssling? and is presented because Fréssling’s 
numerical solution is to be used in the present paper, 


owing to its greater accuracy. 


(3) THE ENERGY EQUATION 


For liquids, the energy equation is (page 606, refer- 
ence +3) 
DT 
plc, = JRV*T + @ (S) 
Dt 
in which / is Joule’s value, 7’ is the absolute temper- 
ature, c, is the specific heat at constant volume, & is the 
thermal conductivity, ® is the dissipation function, DT = 
Dt is the substantial derivative of 7, and V°7" is the La- 
placian of 7. For gases, the work done by pressure 
cannot be neglected even if the dilatation of the ve- 
locity vector is small, and the energy equation is of the 


following form (page 607, reference 3): 


DT —) 
(9) 
Dt 


pit, = JkV*T + ( + 
in which c, is the specific heat at constant pressure. 

If a is defined to be k, pc, c standing for c, for liquids 
and for c, for gases, and if the nonhomogeneous terms 
are disregarded, Eqs. (8) and (9) can be represented by 
the single equation 


DT/Dt = aV?T (10) 


A solution of Eq. (10) is a complementary solution of 
Eq. (8) or (9). The complete solution is obtained by 
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Graphical representation of Eq. (20) for various Prandtl Numbers 


adding to it a particular solution. In the following, 
the complementary solution will be obtained first, and 


the particular solution will be discussed separately. 


(4) DIFFUSION FROM A POINT SOURCE OF HEAT 


Neglecting the temperature change due to pressure 
and viscous dissipation, the temperature rise in the 
fluid due to a point source of heat situated at the stag- 
nation point can be calculated from Eq. (10), which for 
steady flow and for the particular coordinates chosen 
assumes the form 

or or 
d w - 
or Oz 


1 


oT Lol OT 
a a ae = pee (11 
or r or Os? j 
For small thermal diffusivity the last equation can be 
written for the thermal boundary layer as 


u(Ol’/or) + w(OT'/dz) = a(0*T\ Oz”) (12 
If the ambient temperature is designated as 7%, the 
dimensionless quantity 
6 = (T — To) To (13 
can be conveniently used instead of 7 in Eq. (12): 
u(O6/Or) + w(08/0z) = a(0*6 Oz* (14 


If the boundary is impervious to heat, 08/0 vanishes 
atz = 0. Then multiplying Eq. (14) by 27r and inte 


grating with respect to z, one has 


i O60 , O06 O86 
Qrru dz + Qarw dz = a = 0 
0 or 0 Oz Oz_Jo 


(15) 


since 00/0z vanishes also at z = ~. By partial inte- 
gration and from Eq. (3), the second integral above 


can be written 


; O60 
2 mrrw dz = 2nrw6 + 
0 Oz 0 


oO 
ng, 


nd 


he 
ig- 
or 


be 


on 
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since w = QO at z = Oand 6 = Oatz = Thus Eq. 


15) can be written 


2 ‘ 
2rrué dz = 0 
or 0 


. 


/ 2rrub dz = H (16 


which is an integral condition of thermal continuity, so 
that // is actually a measure of the strength of the point 


or 


source. 
If one puts 


6 = (H/V 28v rr*)g(n) (a7 
Eq. (14) becomes, from Eq. (5), 
—a(f’g + fe’) = g” (18) 


in which o is the Prandtl Number v/a, and Eq. (16) be- 


comes 


| f'g dy = 1 (19) 
0 


Phe solution of Eq. (18) with the boundary conditions 
g(~) = 0, g’(0) = Ois 


> 


g=KkK exp (—o f dn) (20) 


in which A is determined from Eq. (19): 


. 


K-11 = J / exp (~o / fan) a (21) 
0 J 0 


Since f is already given by Fréssling,” the last two equa- 
tions constitute the solution for Eq. (18), and Eqs. (13) 
and (17) give the temperature distribution. The 
function g(n),/A is plotted in Fig. 2 for various values 
of the Prandtl Number. The quality A computed 
from Eq. (21) is plotted in Fig. 3 against the Prandtl 
Number. The dependence of A on o can be expressed 


approximately by the formula 
K = 0.8860°°* (22 
almost up to ¢ = 1,024, the highest Prandtl Number 
considered. 
(5) DIFFUSION FROM BOUNDARY 


If 6 varies as an arbitrary power of r on the boundary, 
one can write 

6 = a,(B/v)”" rg, (n) (23) 

in which 7 is any arbitrary positive number and a, an 

arbitrarily assigned constant. Substitution of Eqs. 


(5) and (23) into Eq. (14) results in the equation 


gn” + offen’ — (n/2)f’g,] = 0 (24) 


Sn 


to be solved with the boundary conditions g,(0) = I, 


ARO x 





O: 
06 | 10 100 1000 


o 


Fic. 3. Graphical representation of Eqs. (21), (22), (26), and 
97 


g,(©) = O, the function / being already given by 
Froéssling.” 

For constant temperature at the boundary, n 0, 
and the condition go9(0) = 1 amounts to taking 


a (7) T/T 
in which 7) is the ambient temperature and 7) is the 
temperature at the boundary. The solution of Eq 


(24) for n = 0 satisfying the boundary conditions is 


. 


gZo(n) = K’ / é "So "dn 25 


See ee | e an 26 
u 


« 


in which 


The function gy given by Eq. (25) is plotted in Fig. 4 
for various values of the Prandtl Number. Since the 


rate of heat transfer from the plate is given by 


b= BS | &(~) Ty — To)g0'(0 
ieee Oz . 0 j 4 : 


2 (~) T; — T)K’ 


it seems desirable to express the quantity A’ as a fune 
tion of the Prandtl Number. The values of A’ cal 
culated from Eq. (26) for different Prandtl Numbers 
are plotted in Fig. 3, from which it is seen that the de 


pendence of A’ on o can be expressed by the formula 
K’ = 0.5390°-* (27 


up to 0 = 258, with slight deviation for higher Prandtl 
Numbers. It should be noted that, for » = 0, the left 
side of Eq. (23) is independent of 7, and Eq. (11) coin 
cides with Eq. (14), the solutions of which are thus 
exact solutions of the homogeneous energy equation. 
Noteworthy also is the fact that the solution for o 

0.7 represents the first term of the Fréssling expansion* 
for the temperature distribution in the laminar bound 
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Fic. 6. Solution of Eq. (32) for ¢ = 0.7. 


ary layer of a blunt-nosed body with a uniform tem- 
perature on the surface exposed to air. 

For n = 1 and 2, the solutions are obtained for the 
Prandtl Number 0.7 only and are plotted in Fig. 5, to- 
gether with the solution for m = 0 and the same Prandtl 
Number. For the purpose of calculating the heat trans- 
fer from the boundary for” = 


that 


1 and 2, it may be noted 


g,/(0) = 0.582. go’(0) = 0.845 
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Graphical representation of Eq. (25) for various Prandtl Numbers 


Since Eq. (24) is linear in g,, numerical solutions for 
higher values of » can be easily obtained if needed. 
Naturally, 6 can vary as a fractional power of 7, and if 
6 varies as a sum of powers or as a convergent power 
series of 7, its solution can be obtained by superposition. 


(6) TEMPERATURE RISE DUE TO ViIscoUS 
DISSIPATION 


The dissipation function is 


a 2 (%*) 49 (“) + * (S)'+ (= 4 x) | 
ld es “\y ~ Loe dz | Or 


which, from Eq. (6), can be written 


® = (28) 


12u6°f"? + 28% pr2f” 9 


From Eq. (28) it can be seen that, for small viscosity 
and values of r which are not too small, the first term 
on the right is negligible compared with the second 
term, especially in the boundary layer where f and /’ are 
small. This situation fails to hold only at large values 
of z. One should remember, however, that the Homann 
flow is not realistic at large values of z, since in reality 
the boundary is finite in extent and the flow is uniform 
(u = O) at large z. Consequently, in actual flows the 
first term on the right of Eq. (28) is negligible at large 
values of z also, so that it can be neglected entirely 
except for small values of ry. Thus the nonhomo- 
geneous terms to be considered in an approximate solu- 
tion of Eq. (8) can be taken to be the second term on the 
right of Eq. (28), which is the term yu(0u/0z)? usually 
used. 

Since the nonhomogeneous term varies as the square 


of r, it should be included in the energy equation for 


6 given by Eq. (23), with m = 2 and Gy» replacing go 
in order to distinguish between the particular and the 
complementary solutions. Substituting Eqs. (13), 
(23), and 

@ = 28% pr2f”? (29) 


into Eq. (8) or (9), one obtains the equation 





TEMPERATURE DISTRIBU 


Geo" + ofGe’ — of’G. + Df"? = 0 (30) 
in which 
D = Bva/ascJTy = Gv"*p dok IT (31) 


k being the thermal conductivity. If the solution of 


G" + ofG’ — of'G+f/" =0 (32) 
is found, then obviously G; = DG, with D given by 
Eq. (31). The solution of Eq. (32) for o = 0.7 with 
the boundary conditions G’(0) = 0 and G(o) = 0 is 


plotted in Fig. 6 and tabulated in Table | 


Since 
6 = a(B/v)r*G. = (B*ur?/RIT))G (33 
and G(0) = 0.726, the temperature rise at the boundary 


(assumed insulated) is 
T. — Ty = 0.72682ur?/kJ (34) 
in which 7°, can be called the eigen-temperature at the 
boundary. 
To determine the thermal effects of pressure change, 


it is necessary to calculate 


Dp/ Dt = u(Op/ or) + u(Op/dz) 


Since 
—(1/p) (Op, Or) = ,(Om,/Or) = B*r 
and, from Eq. (2 


Ow 1 Op “qe 


os en - 
Oz p Oz Oz? 
w being independent of r from the second of Eqs. (5), 


it follows from Eqs. (5) that 
Dp/Dt = 4pB?(f?f" + ff") — pB rf’ 
so that 


ce 9) — (Dp Dt) = $u8?(f?f" A. ff” + af’ 2) 4 
pB*r?(2f" ———_ f 


It can be seen that even if, for vanishingly small vis- 
cosity, one neglects the three terms containing u ex- 
plicitly, the term — p6*r?/’ still remains to be considered, 
as well as that given in Eq. (29). For small values of », 
the value of /’ is small compared with 2/"*. But at 
the edge of the boundary layer the reverse is true. To 
account for the temperature change due to pressure 
variation as well as dissipation, it is therefore necessary 
(and sufficient for small viscosity) to take 


pB°r? (or oc aa F 


as the nonhomogeneous terms. The resulting equation 


corresponding to Eq. (30) is 
70" + afG’ — af’G» + D{f" ? iam (f’ 2)] = () 


Since f’(~) = 1 and f"(~) = O, this equation is in- 


compatible with the boundary condition G.(2) = 0. 
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Fic. 7. Dimensionless plot of a typical pattern of isotherms 








with temperature rise due to viscous shear included 7) = ambi- 
ent temperature, AT = 7 — 7 

4 ~ —_—_—— 

3 

2 

if 

| el od 

fe) 1 

Le) | 2 3 

r 


Fic. 8. Dimensionless plot of typical temperature distribution 
at various radial distances. 7) = ambient temperature, 7, = 
temperature at boundary, 7; — 7) = 0.27 


Thus a realistic solution for the temperature distribu 
tion in Homann flow cannot be obtained for gases if the 
thermal effect of pressure change is to be included. In 
the following sections, the discussion and the conclu 
sions are subject to the limitation that this effect is 


neglected. For liquids, the thermal effect of pressure 
change is exceedingly small, and the neglecting of this 


effect naturally does not impose a real limitation. 


(7) DIscuSSION 


The temperature rise due to viscous shear should be 
superposed on the temperature distribution obtained 
from the homogeneous energy equation. If the tem 
perature on the boundary varies with r”, then for n 
less than 2 the heat transfer is from plate to fluid for 
small values of r and from fluid to plate for large values 


of r. The reverse is true if m is greater than 2. For 
n = 2, the heat transfer is either always from plate to 


fluid or always from fluid to plate, or, if the temper 
ature on the boundary is exactly the eigen-temperature, 
there is no heat transfer between boundary and fluid. 
For constant temperature at the boundary, the tem 
perature distribution is obtained by superposition, 
7 aes ry 7, _ T * 1G 0.796 ° 
ro - a £o() n-G(n) — V.¢2071"g2(n) 
To T\ 
in which 


r;?> = Bur?/RIT 
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TABLE | 
Solution of Eq. (32) for Prandtl Number 0.7 


n G n G 
0.0 0.7263 3.0 0.0263 
0.3 0.7091 3.2 0.0173 
0.4 0.6663 3.4 0.0111 
0.6 0.6058 3.6 0.0070 
0.8 0.5349 3.8 0.0043 
1.0 0.4595 1.0 0.0025 
3 0.3846 we 0.0015 
1.4 0.3139 4.4 0.0008 
1.6 0.2501 1.6 0.0004 
1.8 0.1945 4.8 0.0002 
2.0 0.1478 5.0 0.0001 
2.2 0.1097 5.3 0.0001 
2.4 0.0795 5.4 0.0000 
2.6 0.0563 5.6 0.0000 
9 @ 


0.0389 


and the last term involving gs is added to achieve unt- 
form temperature on the boundary. Writing A7’ = 


T — To, ATo = 7; — To, one has 


AT/T) = Ago(n) + n7[G(n) — 0.726g0(n) | 


from which the isotherms can be drawn in the plane of 


In Fig. 7 the pattern of isotherms corre- 
The temperature dis- 


r, and 7. 
sponding to A = 0.2 is shown. 
tribution at various values of r; is shown in Fig. 8, in 
which it is seen that the temperature rise due to dissi- 
pation becomes more and more pronounced as 7” in- 


creases and is predominant for large values of 7. 


(S) CONCLUSIONS 


In laminar stagnation flow with axisymmetry, the 
temperature distribution in the fluid due to a point 
source of heat situated at the stagnation point is given 
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by Eqs. (13), (17), (20), and (21). The function g(7) 


is plotted in Fig. 2. For a general class of concentric 
heating of the boundary, the temperature distribution 
is given by Eqs. (13) and (23) and by the solutions of 
Eq. (24) for various values of 7. The solution of Eq. 
(24) for constant temperature on the boundary (7 = 0 
is given by Eq. (25) and is shown graphically in Fig. 4 


The de 


pendence of the constants A and A’ in Eqs. (20) and 


for various values of the Prandtl Number. 
(25) on the Prandtl Number is shown in Fig. 4, the 
functional relationships being described approximately 
by Eqs. (22) and (27). Forn | and 2, the solutions 
of Eq. (24) are given graphically in Fig. 5 for Prandtl 
Number 0.7. Temperature rise due to dissipation is 
given by Eqs. (13) and (33), in which the function G is 
given graphically in Fig. 6 for Prandtl Number 0.7, 
and is tabulated in Table 1. This temperature rise is 
to be superposed on that which results from the solu- 
tion of the homogeneous energy equation. With the 
thermal effect of the change of pressure neglected, a 
typical pattern of isotherms for constant temperature 
at the plate is shown in Fig. 7, and the vertical tem 
perature distribution at various distances from the 


stagnation point is shown in Fig. S. 
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SUMMARY 


An approximate load-axial deformation relationship is derived 
for pin end columns well into the inelastic range. A series of 
tests for L/?t from 34 to 159 are shown for correlation with this 
ipproximation 

The concept of ‘‘plastic-dynamic buckling”’ is introduced to 
emphasize the strain-rate energy dissipation in early plastic 
buckling of columns This physical concept is applied to a 
generalized structure whose maximum load-carrying capacity is 


reached in the inelastic range 


SYMBOLS 


| cross-sectional area of a member 

Kk = modulus of elasticity 

h depth of a cross section of a member (for rectangular 
section 

I moment of inertia 

1 radius of gyration 

8 = length of a member 

Ler = the critical length of a column (transition length 

2s the axial load on a member 

P-, = the critical load on a column ( Euler’s load 

») stress 

Si yield stress (lower yield stress 

ds = anincrement of length along the longitudinal axis of the 


column 

x the coordinate length taken parallel to the length of a 
member 

y = the deflection coordinate perpendicular to the length of a 
member 

S the coordinate length oll a cress section measured from 
and perpendicular to the axis of minimum moment of 
inertia of a cross section of a member 

A, = change in length of a member 

Ar = the deflection at the center of a column 

the distance from the neutral axis to the yielding fiber of 


a cross section 


INTRODUCTION 


” MUCH OF THE LITERATURE of strength analysis a 
great deal of emphasis has been placed on the maxt- 
mum axial load a column will withstand. This load 
will usually be the controlling factor in a statically 
determinant structure where failure of one member 
means failure of the whole structure. However, in a 
redundant structure, a column may or may not be at 
its maximum load at the instant the maximum load for 


Received July 16, 1953 

* The author wishes to acknowledge the help and guidance of 
Prof. J]. A. Van den Broek, the encouraging efforts of Prof. T. A 
Hunter and Mr. John Rowley, and the patient assistance of his 
wife 

+ Formerly with Department of Engineering Mechanics of the 


University of Michigan 


PARIS 


Lehigh University 


the entire structure is reached. The axial deformation 
of the column may be bevend the point of maximum 
load or may not have reached it at all on its load 
deformation curve. 

In order that an analysis of structures might be 
formulated in which this fact can be taken into ac 
count, it is necessary to derive the load-deformation 
relationship for each member of the structure. From 
these individual member relationships an overall load 
deformation relationship for the entire structure can 
be found. 
all relationship will be the true criterion of the strength 


The maximum value indicated by this over 


of the entire structure. 

The object of this paper is to determine an approxi- 
mate load-axial deformation relationship for pin-jointed 
columns as it might be used in an analysis of pin 
jointed structures. The load-axial deformation rela 
tionship for tension members is now known as the 
stress-strain curve of the material. The addition of 
the load-axial deformation relationship for compression 
members will complete the requirements for solving 
other similar structural 


truss problems and many 


problems. 


THEORETICAL ANALYSIS 


The load-axial deformation relationships for com- 
pression members in the elastic range were obtained 
by Hooke and Euler over 200 years ago. If the column 
is shorter than critical length,{ it will load linearly with 
no buckling until yield stress, Sj, is reached, whereupon 
it will crush and or buckle plastically (Fig. la). A 
column longer than critical length will load linearly 


t Critical length is the length at which the Euler load (P? 
r?EI/L?) and the crushing load (P SA 


eV EI/S,A 


ire equal, or / 





__ P#S.A_ Ded 














Aa Aa 
(a) L<Ler (b) L>Ler 


Fic. 1. Elastic theory for load deformation curves 
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Fic. 4. Assumed plastic stress distribution. 


with no buckling until the load approaches Euler's 
load. On approaching Euler’s load, the column will 
buckle elastically with no appreciable change in load 
as axial deformation increases (Fig. 1b). As the first 
fiber of the member yields, the buckling is said to be 
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plastic. This transition phenomenon is accompanied 
by a decrease in load as axial deformation increases. 
The linear elastic portion for all lengths is expressed 
by Hooke’s Law as 
A, = PL/AE 1) 
Euler’s load is 
2 —_ 2h 2 ‘ 
P, = gl /L (2) 


These relationships are shown in Fig. 1. 


LoOAD-AXIAL DEFORMATION 


RELATIONSHIP 


THE PLASTIC 


After yielding, the plastic column will require only a 
single relationship to represent the system, since only 
one continuous physical change (plastic buckling) is 
occurring with increasing axial deformation. 

The axial deformation of a column in the plastic 
range may be represented as the sum of two causes 
the shortening due to direct stress and that due to 


flexure. This may be expressed by (see Fig. 2) 


| B 7 
A, = ie J (ds — dx) (a) 
AE we 


The integral may be expressed as 


” l ™ dy 
(ds — dx) = dx (4) 
J0 2 Jo dx 


On substituting this into Eq. (5), 


PL | 


“ih 2 
4 | (o) j 
= dx 
AE 2/0 dx 
In order that Eq. (5) may be integrated, it is neces- 


sary to know the deflection curve of the column. Lin! 
has shown that, for early plastic buckling (when only a 


Ay 


small portion of the cross section is in the plastic range), 
the deflection curve may be closely approximated by a 
sine wave. Also worth noting is that this integration 
is of the Rayleigh type 
approximate slope form yields a good approximation 


that is, integration of an 


as a result. 
Considering these facts, let the elastic curve be as- 


sumed sinusoidal. Then 
y = Arsin (rx/L) (6) 
Differentiating, 
dy/dx = (Arm ’L) cos (rx/L) (7) 


Substituting this into Eq. (5) and integrating, 
Aa = (PL Ak) + (a A?7? 1L) (S) 


In order for Eq. (8) to be revised so that it will be a 
function of P only, an expression must be found for 
Ar as a function of P. This may be done by consider- 
ing the stress distribution at the central cross section of 


the column (see Fig. 3). From statics 





ied 


sed 





LinNi?T 


P= SSdz; M = P (Ar) = S'Szdz (9) 
If the type of distribution as shown in Fig. 3 is to be 
assumed, p must be evaluated if Eqs. (9) are to yield 
Ar as a function of P. The expression for p must in 
clude curvature that would require a more exact de 
flection curve than the sine approximation. A more 
exact deflection curve may be found by the method of 
Lin,' which could be differentiated to obtain the curva 
ture. The complexities involved would certainly be 
fruitless, since the net result would be A, as a compli 
cated function of P. Since this relationship will be used 
in further algebraic processes in the attempt to solve 
individual problems, a simple function is essential. 
It seems that a rational approach to this problem lies 
within the works of Van den Broek® in the assumption 
of a rectangular stress distribution (Fig. +). Combin 
ing Eqs. (9) and solving for Az, 


Ar = S Sz dz SS dz (10) 


Substituting Eq. (10) into Eq. (S), 


PL 2 "Ss dz\? 
Aa = * (- ) (11) 


AE * iL \ SSdz 


With the assumption of rectangular stress distribution, 
the integral portions of Eq. (11) may be found directly 
in terms of P and the constants of the given column 
cross section. 

For rectangular cross sections 


(12) 


A, = " 
Pp SA 


PL wr? [h (S;A P ) 
AE . iL L4 ( 

Eq. (12) is the approximation for the load-axial de 

formation relationship for a plastic column of rec- 

tangular section. It may be applied to many other 

sections such as I-beams and angles by calculating an 

equivalent rectangle for the section, provided that 


local buckling does not take place. 


EXPERIMENTAL ANALYSIS 


In the derivation of the plastic load-axial deforma- 
tion relationship a great many assumptions have been 
made. It is therefore necessary to examine these 
relationships experimentally. This presents the prob- 
lem of devising an accurate testing method that will be 
sufficiently accurate in the range of large deformations 
encountered in the presented theory. Time is a most 
important element during such tests, since strain rates 
affect the yield stress of materials. This time factor 
requires simplicity in each test. 

One-half-inch-square, annealed, hot-rolled steel was 
selected as the material to be tested. All specimens 
were cut approximately to length and annealed in the 
same heat. Five specimens were tested in tension and 
five in compression to determine the stress-strain rela- 


tionship of the material. An average yield stress of 
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Fic. 5. Cylindrical pinned joints 
approximately 35,000 Ibs. per sq.in. was found. The 
greatest deviation in yield stress among the ten tests 
was 3 per cent. No appreciable change in the yield 
stress was recorded up to 0.02 in. per in. of strain. The 
material exhibited the ideal plastic properties in much 
the same way as ordinary structural steel, the only 
appreciable difference being that residual stresses were 
eliminated by annealing. 

A most important fact is that no effort whatsoever 
was made to obtain straight specimens. To this fact 
alone might be attributed the deviations in attaining 
the crushing load (P = 5S,A) and the Euler load (P = 
r°EI/L*) in the following test Certainly, 
initially crooked columns are more closely related to 


results. 


those in actual structures than are straight ones. 
Because of this fact the rationality of the critical load, 
Euler's load, has been a topic of considerable con 
troversy. Initially were herein 
used to show that, with columns in which the elastic 
analysis may be in error due to this crookedness, 


crooked specimens 


the presented plastic theory will rapidly converge to 
agreement with an idealized theory. This is due to the 
fact that in the plastic range the order of magnitude of 
deformations is much larger than the order of magnitude 
of the initial crookedness. One might speculate further 
and say that, therefore, a more rational general column 
theory might result from a large deformation theory of 
columns. 

The next consideration in the tests was in attaining a 
perfect pin-ended column. By placing the specimen 
in a slot machined precisely in the center of a half 
cylinder and letting the cylinder roll freely on ground 
surface plates, a pin end was attained (see Fig. 5). 
The line of action of axial load will always be through 
the line that is the centerline of the cylinder. The 
direction of the buckling is then controlled so that 
transverse deflection gages may be employed with little 
difficulty. The accuracy of this technique may best be 
illustrated in the attainment of Euler loads and crushing 
loads in the author’s previous work.‘ Fig. 6 shows the 
results of these tests in which the same cylindrical 
heads were employed. 

The surface plates were mounted in the testing ma 
chine and shimmed parallel with an accuracy of 0.0001 


in.in Sin. of length. The testing machine chosen was a 
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ends of the specimen. A O.0001-in. dial gage was 
PREVIOUS TESTS employed to measure directly this change in length. 
wine CYLINDRICAL Two 0.001-in. dial gages were used to measure the 
: a Se er oe transverse deflection of the center of the specimen 
. P=SA = BARR STOCK ge oA soar 
the Fig. 7 shows the specimen, heads, and axial deflection 
a gage. Fig. S shows the apparatus in the machine ready 
= 
pd to run a test. 
a 
. P = Per DISCUSSION OF RESULTS 
In Figs. 9 to 12 the curves represent Eq. (i2) pre- 
sented in the section on ‘Theoretical Analysis."’ The 
fe) — : points marked (0) are the average values of the test 
0 50 100 Li 150 data. The maximum deviation in the test data in the 
plastic range was less than 7 per cent for all but the 


Fic. 6. Previous results with cylindrical heads for column tests. 





Fic. 7. The specimen, heads, and axial deflection gage. 





Fic. 8 


Setup in the testing machine. 


Riehle screw-type machine with electronic controls. 
This made possible accurate control of the rate of axial 
deformation on which strain rates are dependent. 

To complete the experimental apparatus, an axial 
deflection gage the axis of the 
“pinned joint’? was on a line (the centerline of the 
cylindrical heads) which passed through the specimen, 
the specimens could be gage-punched at these points. 
A gage much like the ordinary mechanical strain 
indicator was built to measure the change in distance 
between these gage-punched “pin joints’’ at opposite 


was devised. Since 


34.6. 
deviation in 


specimens of L/t = This compares favorably 


with a maximum attainment of Euler 


loads of approximately 10 per cent. 
less than 7 per cent deviation of test data in the plastic 


The exception to 


range was short specimens of 1/7 = 34.6 (shown on 
Fig. 9). 
and continued to 
crushing without buckling is a direct violation of the 
No way of predicting this 


These specimens reached the crushing load 
deform without buckling. This 
assumptions of the theory. 
phenomenon, either to whether it would take place or 
for how long a specimen would continue to crush with- 
out buckling, was found. At for the short- 
column range, the theory presented may be considered 


least 


as the minimum load condition of a load-axial deforma- 
tion relationship for columns of Li less than 50. Ina 
structure, any errors due to use of this theory will be on 
the safe side. 

For the longer column specimens of L 7 of 76.5, 118, 
and 159, the correlation between theory and test is 
somewhat astounding. The idealized plastic theory, 
with all its assumptions represents the true physical 
picture with more accuracy than does the idealized 
elastic theory. Not only does it do this, but the plastic 
theory becomes more accurate as deformations in- 
crease. Figs. 10 to 12 are evidence of these results. 
Certainly the desired plastic load-axial deformation 


relationship is well approximated by the theory as 





presented. 
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Fic. 9. Load-deformation curve for L/i = 34.6 
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L/i= 763 
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Fic. 10. Load-deformation curve for 1/1 = 76.3. 
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Fic. 11. Load-deformation curve for L/z = 118. 


CONSIDERATIONS IN COMPARING THEORY 


AND TEST 


FURTHER 


It should be emphasized that the derived theory 
is a statical analysis. No considerations have been 
made for strain rates that must occur in testing as 
well as in actual structures. In Fig. 13, which is an 
enlarged plot of the first part of Fig. 10, it can be seen 
that the experimental data (dashed lines) lie well above 
the statical theory in the early plastic buckling range. 
This plot is for the specimens of L/i = 76.3, the spect- 
mens that are closest in length to the critical length. 
The differences between the static theory and data 
recorded may be attributed to the effect of strain rates 
in dynamic buckling. this dynamic buck- 
ling occurs in the first portion of the plastic range, 
it has been given the name ‘‘plastic dynamic buck- 


Since 


ling.” 

Fig. 13, which best shows ‘“‘plastic dynamic buckling, ”’ 
is a plot of load vs. axial deflection with deflection 
measured at the point of the load and in the direction 
of the load. Therefore, the areas on Fig. 13 represent 
A concept of plastic dynamic buckling in 


energy. 
terms of the dissipated energy may be arrived at by 
careful examination of the areas of this graph and 


physical considerations of these tests. 
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THE CONCEPT OF PLASTIC DYNAMIC BUCKLING 


The testing machine itself must be considered in 
the concept. The mechanism of a screw machine is 
such that A, will always be increasing unless the re- 
verse button of the machine is pushed. Since, in the 
tests, the machine was set to a slow rate and left on 
throughout the test, A, was constantly increasing. 
This was true in all observed data. Another consider- 
ation to be made with respect to the testing machine 
is the slope of line AD in Fig. 13. The time required 
for the specimen to proceed physically from A to D 
is extremely small; thus little additional axial deforma 
tion may be attributed to the turning of the screws of 
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Fic. 14. A P vs. L/i plot showing the limit of dynamic insta 
bility. 


the testing machine. The comparatively large addi- 
tional deformation recorded from A to D is due to the 
release of strains within the testing machine itself. 
If the testing machine is considered to act as a linear 
spring, the line AD will be a straight line. With these 
concepts in mind, a straight line has been drawn from 
A to D and considered actual test data. Evidence 
that this linear spring effect for the testing machine is 
correct may be found from examination of the slope of 
line AD for each specimen tested, which reveals in all 
cases that these lines have nearly the same _ slope. 
Since this effect was observed for cases in which the 
loads at points A and D were of different magnitudes, 
the close similarity of slopes indicates that the testing 
machine did act as a linear spring. 

The statical theory indicates that, as plastic buckling 
occurs, the load-axial deformation relationship should 
follow the curve ABC. Since, in a testing machine, as 
in a structure, the axial deformations are constantly 
increasing, this curve cannot be followed in an actual 
column. The energy in the area ABCD must be dis- 
sipated by some means other than those considered 
in the derivation of the “‘statical’”’ theory. 

There are two possible means by which this energy 
may have been dissipated. They are: (a) the plastic 
strain rates involved in sudden buckling that would 
increase the yield stress of the material during sudden 
buckling, and (b) any energy that might be dissipated 
by accelerations of the mass of the column. By obser- 
vation of many tests of columns and material properties, 
it becomes apparent that the energy dissipated by 
acceleration is small. The area ABCD on Fig. 13 
has therefore been labeled as strain-rate energy dissi- 
pated in plastic dynamic buckling. The conclusion to 
be drawn from this phenomenon as shown is that the 
plastic dynamic buckling of a column is a process 
directly dependent on strain rate in the plastic portion 
of the column. The reason for suddenness of plastic 
buckling is that the statical column cannot store or 
dissipate the amount of energy which has been put into 
the column during the elastic loading process. Specula- 
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tion might well be made on the thermodynamic effects 
of the nonreversible process, plastic buckling of columns. 
Consideration of such effects is beyond the scope of 
intention of the author. It is introduced as a point 
to be considered in further investigations. 

Another interesting observation from Fig. 13 is that 
the straight line AB would represent the line AD for a 
column tested in an infinitely rigid testing machine. 
The load corresponding to the load at B is plotted 
versus length (or L/z, since 7 is constant) in Fig. 14. 
The experimental loads corresponding to points D for 
for each specimen are shown marked (0). From con- 
sideration of Fig. 13, the line on Fig. 14, which is deter- 
mined by point B, is the maximum load at which dy- 
namic buckling could cease. The experimental points 
support this observation rather well. 

Although the load-axial deformation relationship 
derived has treated the statical cases only, it also sheds 
light on the mechanism of dynamic buckling. Enough 
so that it is possible to consider the transition of the 
discussion from a plastic column in a testing machine 
to that of a column in a structure. 


APPLICATION TO THE GENERALIZED STRUCTURE 


Let it be assumed that a pin-jointed column is a 
member of a general redundant structure. Let the 
structure be of such a nature that it is completely 
solvable by present methods, except for the action of 
the column as a part of this structure. For simplicity 
in this illustration, assume that the structure is loaded 
with a concentrated force along the centerline of the 
unsolvable column. Let the column be removed from 
the structure and replaced by its reaction P (see Fig. 
15). The deflection of the structure without the 
column can now be solved in the form 


A= Ai(Q — P) 


The nature of this relationship is a load-deformation 
relationship of the structure. Since A and A, (as 
applied to a column) are equal, the load-deformation 
relationship for the column may be derived from the 
theory presented in this paper as 
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Fic. 15. A generalized structure. 
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LIMIT DESIGN 


A = fo( P) 


Che sum of these relationships for a given A = A, will 
vield the load that the structure will sustain at the 
deflection A). 

The simplest method to solve for the maximum load 
the structure will sustain is shown graphically in Fig. 16. 
Fig. 16b shows a typical load deformation curve for a 
column near critical length, the worst possible type of 
column for this illustration, since, compared with other 
lengths, its load drop in the plastic range is more pro- 
nounced. Even for a structure with this near critical 
length column, the maximum load carrying capacity 
may actually come after the column has entered the 
plastic range. (This is clearly shown in Fig. 16c.) 
The areas of strain-rate energy may be found in the 
same manner as was shown for a column in a testing 
machine, but a simpler method may be reached by 
inspection of Fig. l6c. It is usually the nature of the 
load Q to be a monotonically increasing load; there 
fore the dashed line rather than the solid line must be 
followed in the actual load deformation path. This 
explains the sudden buckling as noted on this illus- 
tration. Thus the maximum load for the structure of 
this illustration is attained after plastic dynamic buck- 


ling has occurred in the column. 


CONCLUSION 


Although many structures may not exhibit this 
maximum with a column in the plastic range, certainly 
others will. Therefore, the limit design of columns 
in structures must be considered if the structural 
analysis is to be complete. 

Recently, a great deal of emphasis has been placed 
on blast analysis of structures. A blast may cause 
sudden large deformations in structures. The energy 
of the blast dissipated in the structure will be in part 
dissipated in the compression members. The theory 
presented in this discussion will afford a means of con 
sidering the large deformations and also for the energy 
dissipation in columns due to such dynamic loading. 

It is worthy to note that equivalent length concepts 
will allow this discussion to be extended to include 
columns with end fixity. Similar treatment may be 
given to beam columns in the manner that pin-jointed 
columns were treated in the ‘theoretical analysis”’ 


of the presented theory. 
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Fic. 16. Column action in the generalized structure. 


In summary, herein is presented a new concept of 
failure of redundant structures containing columns. 
This concept may, in part, explain many of the in 
stances where testing has indicated overdesign. This 
new theory will be applicable in cases where minimum 
weight and strength are a primary consideration, 
although limitations of use are realized due to large 
deformations. 

REFERENCES 

‘Lin, Tung-Hua, J/nelastic Column Buckling, Journal of the 
Aeronautical Sciences, Vol. 17, No. 3, p. 159, March, 1950 

2 Bleich, Friedrich, Buckling Strength of Metal Structures, 
Ist Ed., pp. 34-41; McGraw-Hill Book Company, Inc., 1952 

3 Van den Broek, John A., Ductile Equilibrium Column Theory, 
Overdruk Voordrachten No. 5-1949 

4 Paris, Paul C., Testing of Columns with Uniformly Distributed 
Transverse Loads, The Michigan Technic, Vol. 71, No. 2, p. 14, 


November, 1952 











Design and Performance of an Adjustable 
Two-Dimensional Nozzle with 
Boundary-Layer Correction’ 


G. H. BACKER? 
Applied Physics Laboratory, The Johns Hopkins University 


SUMMARY 


A two-dimensional nozzle with boundary-layer correction to 
both the shaped and unshaped walls was designed with an adjust- 
able throat to operate at the nominal Mach Number 3.106 
and a static jet pressure of 85 microns of Hg (0.0236 Ib. per sq.ft.) 
The nozzle was evaluated in a low-density free-jet wind tunnel 


(No. 2) located at the Low Pressures Project, University of Cali- 
fornia, and met substantially all of the design conditions with the 
exception of the desired free-stream cross section. Experimental 
results verified that it would operate over a limited Mach Num- 
ber range without producing observable disturbances to the uni- 


form exit flow 


NOMENCLATURE 


A = boundary-layer constant 

B = boundary-layer constant 

b = width of ideal nozzle (in. ) 

H half height of shaped walls of an ideal nozzle at exit 
(in. ) 

h = half height of shaped or parallel walls of ideal nozzle 
(in.) 

h* = halfheight at throat (in. ) 

LL = axiallength supersonic ideal nozzle (in. ) 

L, = length of ideal shaped wall (in. ) 

K = coiistant 

M == Mach Number 

Ms; = Mach Number at boundary-layer seam 

p = pressure (microns of Hg) 

pe = test chamber static pressure (microns of Hg) 

pi = impact pressure (microns of Hg) 

pj = static jet pressure (microns of Hg) 

fit = ideal nozzle static pressure (microns of Hg) 

po = stagnation pressure (microns of Hg) 

p2 = surface pressure on cone (microns of Hg) 

u = velocity-axial (ft. per sec.) 

us = velocity-axial-boundary layer seam (ft. per sec.) 

“ur = velocity-axial-ideal at throat (ft. per sec. ) 

WW = weight rate flow of air (lbs. per hour ) 

XY = axial distance from nozzle throat (in. ) 

XY’ = distance along shaped wall from nozzle throat (in.) 

y = vertical distance from nozzle axis (in. ) 

6 = boundary-layer thickness (in. ) 

6* = displacement thickness (in. ) 
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6* = 6*/L (in.) 
n = 9/6 
Mo = viscosity of air at stagnation temperature 
(Ibs. sec. /ft.?) 
p = density of air (Ibs. sec.?/ft.*) 
ps = density of air at boundary layer seam (lbs. sec.?/ft 


INTRODUCTION 


— Low PRESSURES PRoJeEcT of the University of 
California has been engaged over a period of years 
in producing supersonic flows in a low-density gas.'~® 
To accomplish this, different types of supersonic nozzles 
have been tried in an effort to produce a sufficiently 
large uniform flow that would permit accurate aero 
dynamic tests of models such as cones and flat plates. 
Of all of the configurations tried, the most successful 
up to date has been the axially symmetric type." 4 
With these nozzles, a considerable wealth of aerody- 
namic data has been obtained with respect to their 
flow and pressure characteristics. Two serious limita- 
tions to the overall usefulness of these nozzles have been 
uncovered from an analysis of these data. The experi- 
mental results show that small circumferential wall de- 
The 
increments of these waves may be small, but the total 
of their effects when focused at the apex of the cone 
produces a disturbance that considerably alters the 
flow and pressure characteristics of the nozzle, espe- 
cially along the axis where the apexes are located. At- 
tempts have been made to calculate the position of the 
wall deviations for the purpose of reworking the nozzles 
Only lately has this procedure been 


viations produce conical sheets of Mach waves. 


to remove them. 
successful, since the changes are so minute as to be prac- 
tically beyond the limits of the normal fabricating 
Another serious limitation is the absolute 
Once installed in the 


methods. 
fixed geometry of the nozzle. 
wind tunnel, there is no possibility of altering the geom- 
etry, especially while it is in operation, to produce a 
variation in Mach Number. 

A two-dimensional type nozzle should, in theory 
(and to a certain extent in practice), eliminate the two 
limitations just considered for the axially symmetric 
type. In the formative days of trying different types 
of nozzles, a two-dimensional fixed geometry type was 
tried without success.* The major difficulty in this 
case was that no boundary-layer correction was made, 
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and consequently the entire exit was filled with bound- 
ary layer and shock waves. To overcome this diffi- 
culty and to attempt removing the limitations of the 
axially symmetric type, a variable geometry two- 
dimensional nozzle with boundary-layer corrections to 
both pairs of walls was designed and evaluated. 


ANALYSIS 


Since a prerequisite of any supersonic nozzle is a satis- 
factory isentropic core, a two-dimensional (parallel 
sidewalls) ideal core must be selected or calculated. 
Because the emphasis of this nozzle design is on the 
geometry and the viscous corrections, no attempt was 
made to calculate the core, but instead a designed core 
known to have produced satisfactory results under test 
was used. Data for such a nozzle core were obtained 
from the NACA Ames Aeronautical Laboratory. 

In order to secure a subsonic monotonic velocity 
increase in the region prior to the throat of the nozzle, 


ft me St awa 
su" ay — Us 5 1u ay 
dx J, ° . "dx Je - , 


By letting 


an analysis of a type devised by Tsien® was employed 
to secure the desired wall curve for a portion of the sub 
sonic inlet. Fig. 1 shows the two resulting subsonic 
contours used. 

The dimensions of the isentropic core based upon a 
maximum weight flow of 0.75 Ib. per hour at a design 
static jet pressure of 85 microns of Hg are as follows 


(a) Throat Dimensions 


Half height, h* 0.0823 in 

Width, } = () 7645 in 
(b) Exit dimensions 

Half height, 7 = ().3822 in 

Width, b = ().7645 in 
(c) Length (supersonic), L = 2.7433 in 


The cross section of the ideal exit is square, and the 
sidewalls are parallel. 

The problem of predicting the necessary nozzle con 
tour changes to accommodate the viscous rather than 
the ideal flow of a low-density gas is formidable and, 
at best, represents only an approximation. A procedure 
commonly used for this calculation of the viscous effects 
is to assume that the von Karman momentum integral 
boundary-layer equations apply. Using these, the dis 
placement thickness along the shaped and unshaped 
walls is determined, and the sum of these two is applied 
to the shaped walls only, leaving the unshaped walls 
parallel. The additional opening provided by the dis 
placement of the shaped walls will allow the proper mass 
flow to pass, but it will change the configuration of the 
ideal core at the exit. Where nozzles are operated at 
higher densities, the displacement thicknesses are 
sufficiently thin as not to change the exit geometry 
by any significant amount. 

In applying this same method to the low-density 
gases, it is found that the displacement thicknesses can 
not be classified as thin. As has been previously 
stated, the fixed geometry nozzle that was tried failed 
because the viscous flow filled the nozzle. However, 
some estimate as to the additional displacement must 
be made. Schaaf® evolved a method applicable to the 
axially symmetric nozzles which incorporates the von 
Karman momentum integral method modified for 
compressible flow. Since this method is presumed to 
be somewhat insensitive to the type of velocity dis- 
tribution in the boundary layer, a uniform linear vari 
ation is used. It is also assumed that the Prandtl 
Number = | through the boundary layer, which iden- 
tically determines the energy equation as constant 
and equal to free-stream energy. The Goldstein’ form 
of the boundary-layer equation used for this analysis 


is as follows: 


dp < du 
a- — I h dy — poh (I 
dx Jy dy} watt 








n = u/u, = y/6, dn = dy/b, du/dy = u,/6 


then 


dé _ | (d/dx) (psus"hA ) 


6 
dx p.u,"h(B — A) 


Since the displacement thickness (6*) is more useful to 


this analysis, it may be obtained by 
6* = 6(1 — B) (4) 


Solving dé, dx by numerical integration and converting 
to 6* gives the displacement thickness along both the 
shaped and unshaped walls. 

No difficulty arises in applying this solution to the 
axially symmetric case, since the displacement thick- 
In the 
two-dimensional low-density flow nozzle, it is not pos- 


ness is applied uniformly to the circumference. 


sible to follow the practice of displacing the shaped 
walls an additional amount to take care of the boundary 
layer on the parallel walls, since the uniform flow exit 
configuration will be so greatly changed because of the 
thick boundary layers as to become useless for test pur- 


poses. A study of the solutions for the axially sym- 
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“I dp du, 
dn, B= pu/ psu, dn, = — yl; 2 
0 dx dx 


(d/dx) (psusB) (d/dx) (us) bo 
bpsush(B — A) 


_ (3 


u,(B — A) bp.u,(B — A) 


metric boundary layer have shown that the rate of its 
growth with length is nearly a linear function. Since 
this same type of solution is used for the parallel walls, 
it was decided to diverge them uniformly to account 
for the boundary layer. In order to minimize the 
possible Mach waves at the start of this divergence, a 
large radius may be incorporated. The displacement 
thickness was applied as calculated to the shaped walls, 
and then to make sure of a completely faired curvature, 
differences in ordinate were adjusted to the third de- 
rivative to secure a smooth contour. 


DESIGN 


One of the advantages of a two-dimensional type 
nozzle is the utilization of geometry changes while under 
test. In this nozzle a rather unique sliding arrange- 
ment, dictated in part by the required vacuum sealing, 
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of the shaped blocks was incorporated to make possible 
area ratio changes that, neglecting contour changes, 
gave a variable Mach Number nozzle. Fig. 2 shows 
the finished nozzle assembly in which the shaped 
blocks are moved vertically with respect to the side- 
walls. 

To determine the expansion characteristics of the 
nozzle, five pressure taps were drilled into one sidewall 
(Fig. 3) and pressure-sensing thermistors were at- 
tached to determine static pressure variations. It was 
possible by moving both shaped blocks in the same di- 
rection to determine static pressure variations at dif- 
ferent vertical positions of the nozzle. By this trav- 
ersing arrangement, a field of pressures was ascer- 
tained of the expansion flow in the nozzle. 
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A most difficult problem in this type of a variable 
geometry nozzle designed for low-density operation is 
that of providing adequate seals to secure vacuum tight 
ness and at the same time permit relative motion of the 
parts so sealed. In the normal density tunnels inflat 
able rubber seals may be used without difficulty, but at 
low densities the gas used for inflation would itself be 
a possible source of a leak into the system, since under 
vacuum conditions the gas would pass through the 
walls of tube. Solid gasket materials that have good 
out-gassing characteristics are required. To seal the 
blocks to the sidewalls, a continuous type of seal was 
devised as shown in Fig. 3. This proved to be adequate 
when sufficient vacuum grease was applied to sidewalls 
and the seal. Although not comparable to the tight 
ness secured in nonmoving joints, it was satisfac- 


tory. 
EXPERIMENTAL FACILITIES AND PROCEDURES 


The nozzle tests were made in the No. 2 wind tun 
nel of the Low Pressures Project, University of Cali 
fornia.® 

To explore the flow conditions at the exit of the 
nozzle, a traversing mechanism having three degrees 
of motion was installed in the test chamber. A conical 
static tube and a source-shaped impact tube were at 
tached to this mechanism, and pressures were read by a 
precision oil manometer* connected to the system. 
The thermistor pressures were read by means of a re 
cording potentiometer. A mercury McLeod pressure 
gage was used to calibrate the precision oil manometer 
and the thermistors. 

The Mach Number, ./, and the static pressure, ps, 
were determined from the impact and cone traverse 
readings. Using a procedure outlined by Kane,’ it 
may be shown that the Mach Number, .\/, and ratio of 
cone pressure, p2, to static pressure in free stream, Pp, 
are functions of the ratio of measured impact pressure, 
pi, and cone pressure, po. By using tables, the com 
puted ratio, p: pe, may be converted to \V/ and p,. 

The above procedure is based on the assumption of 
ideal compressible flow. Corrections to both the im 
pact and static pressure readings are required to ac- 
count for viscous effects. One set of data (h* 
0.082 in., FR = 0.75 lb. per hour) was corrected for 
these effects, using experimental data from reference 
10 for impact pressures and the predicted changes 
from reference 11 for cone pressures. 

The displacement thickness for the design condh- 
tion was computed from the experimental vertical 
Mach Number distribution over the nozzle exit by 


” VK ) , 
_— l a dv (5 
. | ( 1+ 02M — kK)” 


where 
K = (AM? 4+ 0.2M,2M?)/(M,? + 0.2M2M 
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EXPERIMENTAL RESULTS 


Figs. 4, 5, 6, 7, and 8 and Table | show the experi 
mental results. 
as the mass flow decreases, the viscous region in 
creases until, at 0.55 Ib. per hour, it entirely fills the 
This is not unexpected, since this same result 


The upper curves of Fig. 4 show that, 


nozzle. 
has been common in the axially symmetric nozzles 
where Mach Number variations have resulted from 
changing the mass flow.' The lower curves are of con 
siderable interest, since they show the results of using a 
variable area ratio while maintaining a constant mass 
flow. Two variations not present in an axially sym 
metric nozzle are evident: One is that the boundary 
layer will fill the exit area as the area ratio is increased 
(ratio of exit to throat area) by moving the nozzle 
blocks closer together, while a broad band of uniform 
flow is present as the ratio is decreased. The other is 
that, as expected, the area ratio change causes a vari 
ation in Mach Number. The latter change is further 
illustrated in Fig. 5 and apparently follows the theo 
retical trend though displaced from it. The absence 
of any marked disturbances to the flow as the area ratio 
is changed indicates that it may be used successfully 
over a Mach Number range without resorting to con 


tour changes in the nozzle. It should be emphasized, 
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TABLE 1 
SIDE WALL STATIC PRESSURES 
h = 0.082 inches W = 0.75 lbs/hr 
POSITION MACH Pr P, (Ideal) Pr ( MEASURED) 
5 . saci Po NOZZLE NOZZLE 
Inches’ Inches Microns BLOCKS No } BLOCKS No 2 
Microns Microns 
0.50 0 2.11 .10770 394 382 370 
1.00 0 2.74 .04039 160 155 151 
1.50 0 3.07 .02452 90 103 101 
2.00 0 3.10 .02345 86 87 85 
2.50 0 3.10 .02345 86 7 79 
0.50 0.10 2.14 .10270 376 343 362 
1.00 0.10 2.74 .04039 160 152 152 
1.50 0.10 3.01 .02682 98 101 102 
2.00 0.10 3.10 .02345 86 87 86 
2.50 0.10 3.10 .02345 86 81 79 
0.50 0.20 320 
1.00 0.20 2.68 .04429 162 154 152 
1.50 0.20 2.92 .03071 102 103 104 
2.00 0.20 3.07 .02452 90 88 86 
2.50 0.20 3.10 .02345 86 81 79 





however, that the instrumentation used might not 
sense the smaller flow disturbances because of the size 
of the probes used. Smaller probes might be tried 
to obtain better detail, but the time constants would 
be greater and out-gassing might become a problem 
thus reducing their reliability. 

Fig. 6 shows similar curves for horizontal traverses. 
The changes are of the same nature as for the vertical 
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traverses just considered. Since no change in horizontal 
dimension of the nozzle was possible, it would be ex- 
pected that the boundary layer would not change ap- 
preciably when the area ratio is changed. The Mach 
Number variation, of course, follows the same trend 
as for the vertical traverses. It should be pointed out 
that a serious mechanical limitation of the traverse 
system did not permit extensive horizontal traverses. 
Since the uniform flow region is only about half the 
width predicted, the boundary layer is greater than 
calculated, and a next step would be to alter the side- 
walls of this nozzle using an increased divergence. 
There is, of course, no certainty that this would not 
produce troublesome Mach waves, but from the open- 
ings that were used with the shaped walls this does not 
appear likely. 

The axial traverses showed no definite trend but 
appear to be quite uniform over the range evaluated. 
This trend is indeed different from the serious devi- 
ations experienced in the axially symmetric nozzles.! 
The testing of slender bodies of revolution could be 
accomplished with this nozzle without difficulty. 

The results of the thermistor pressure measurements 
along the sidewall of the nozzle are shown in Table 1. 
Comparison of the measured to the ideal values shows a 
remarkable agreement. When it is possible that these 
pressures may be +10 microns in error, the agreement 


is even more remarkable. Whether these results are 


Boundary-layer profile at exit 


fortuitous or not, it would appear that the low-density 
gas expansion does follow the calculated path with at 
least the same order of regularity as for those used in 
higher density operation. 

The curves shown so far have been based on instru 
ment interpretations assuming an ideal compressible 
fluid. Fig. 7 shows the effects of applying viscous cor- 
rections to the pressure readings. Of greatest interest 
are the upper curves showing the static pressure vari 
ations. By applying viscous corrections to cone pres- 
sure, a modified static pressure curve results which is, 
within experimental error, uniform across the nozzle. 
This follows the boundary-layer premise of no pressure 
change through this region. By re-evaluating the 
Mach Number based on the corrected cone pressure, 
the result gives a centerline Mach Number within ; 
few per cent of the ideal. It should be pointed out that 
the cone correction!! is a prediction based on a linearized 
theory and has as yet no experimental verification other 
than the results shown here. 

The displacement thickness of the boundary layer 
as calculated shows that, for the shaped walls, 6* = 
0.2736 in. at the exit. By experiment, using the Mach 
Number variation for the design condition, 6* = 0.2479 
in. This agreement is good, but it does not tell the 
whole story. Fig. S shows the curves for uu, and 
pu/pu, for both calculation and experiment. It is 
evident that the uniform linear distribution of velocity 
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assumed does not fit the curved profile that results. 
It would appear from the curves that a uniform mo- 
mentum curve would more nearly fit the experimental 
results. The uniform flow, free from the boundary 
layer assumed, is not entirely achieved. 

A comparison of results for the two nozzle blocks 
shown in Fig. | indicated that, within the expected 
experimental error, there is no difference between the 
two Evidently the has 
little effect on the resulting flow of these two configu 


nozzles. subsonic entrance 


rations. 
CONCLUSIONS 


A two-dimensional adjustable area-ratio nozzle with 
boundary-layer correction has been designed to pro- 
duce wave-free uniform parallel supersonic flow at its 
exit. The experimental evaluation of this nozzle shows 
that, although the desired results were achieved only to 
a limited degree, they do indicate that the concepts 
are sound and should be investigated further. 

The divergence of the sidewalls, although not suffi- 
cient to take care of the boundary layer to the extent 
desired, does give a small cross section of practically 
uniform flow. Since this flow is apparently free from 
wave disturbances, a further divergence should be 
investigated to determine just how far it is necessary to 
diverge to attain the desired region of uniform flow. 

The vertical sliding of the shaped nozzle blocks pro- 
duces consequent area ratio changes that result in 
a variable Mach Number nozzle over the range tested. 
Fig. 5 indicates that the experimental results follow 
the theoretical area ratio trends. If pressure instru 
mentation had been corrected for each ratio, it is prob 
able that the agreement would have been even better. 
Of equal importance is the apparent absence of flow 
disturbances for these Mach Number variations. Since 
the nozzle was designed for a discrete Mach Number, 
M = 


istics would show waves appearing for any other Mach 


3.106, two-dimensional isentropic flow character 
Number. This absence may be due to the nature of 
the pressure instrumentation or the manner in which 
it was used. It might also be that boundary-layer 
changes are of such a nature that a true isentropic core 
results. It might also be a combination of the above 
reasons plus others which mask out any evidence of 


flow changes. It would be expected that a larger ver- 
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sion of this nozzle would determine whether a truly 
variable Mach Number nozzle with adequate flow 
characteristics has been achieved. 

A re-examination of the boundary-layer calculation 
using a diiferent concept for the predicted velocity 
variation might give a better estimate of the viscous 


layer. Owen‘ has tried this for an axially symmetric 
nozzle. In this case, a sine variation was used for 
u = f(y). The suggested use of the uniform momen 


tum gradient might also be tried. 
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nen itself responsible for the opinions expressed by the correspondents Contributions should no 
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niet On ‘‘Supersonic Flow About Slender Bodies of REFERENCES 
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t 30, Elliptic Cross Section”’ Kahane, A., and Solarski, A., Supersonic Flow About Slender B 
Elliptic Cre Section, Journal of the Aeronautical Sciences, Vol. 20, No. 8 
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I given in the paper on this subject by Kahane and Solarski! + 
Mach were previously obtained by Fraenkel,? whose paper, however, 
ty ol has had only a limited circulation to date ‘ 
23 In spite of the similarity of their titles, these two papers do not Gap Effect on Slender Wing-Body Interference 
cover precisely the same ground. Reference 1 permits axial ; 
Vind variation of the eccentricity of the cross sections but requires the Duane W. Dugan - 
». 2, body to be pointed and of continuous profile slope, whereas Ames Aeronautical Laboratory, NACA, Moffett Field, Calif 
reference 2 requires constant eccentricity of the cross sections but October 5, 1953 
nary permits the profile slope {and therefore the source strength func 
jects tion S’(s)] to be discontinuous. (Open nose bodies are thus N A RECENT CONTRIBUTION to the Readers’ Forum,! H. Mirels 
included in the treatment of reference 2 Nevertheless, the two pointed out that when the gap effect is included, the solution 
ICs, papers have much in common. Fraenkel applied his general to the lift effectiveness of a slender wing hinged to a cylindrical 
theory to elliptic cones, and he also extended some previous calcu body is unobtamable without recourse to numerical methods 
tres lations* of the wave drag of slender bodies of revolution having It may be interesting to note that the writer has obtained in 
75 straight and parabolic profiles (and pointed or truncated ends closed form the solution to the problem of the loss in lift caused 
23, to bodies of elliptic cross section by a clearance gap between the wing panels and body of a slender 
Simple general expressions for the lift, induced drag, and pitch wing-body combination having the body at zero angle of attack 
per ing moment of elliptic bodies (with constant eccentricity and dis and the wing deflected In comparing the results with the effect 
ects continuous profile slope) are also given in reference 2. In this Of 4 clearance gap upon the lift of a slender wing-body combina 
the connection it should be mentioned that the remark on the inci tion having wing and body at the same angle of attack, it is not 
PP dence problem made in reference 1 (on page 521 namely, that worthy that the loss of lift, percentage-wise, is considerably 
smaller than in the latter case In either case, an increase in 
pre Ap/q = (Ap/q)a=0 + (AP/q)a body-diameter-span ratio results in a greater loss in lift effective 
‘ali is dangerous, because of the fact that the pressure coefficient has — ness for a given gap-span ratio 
23 nonlinear terms and may lead the reader to the well-known error In the limiting case of zero gap between deflected wing panels 
0) of omitting the coupling terms arising from the thickness and in and noninclined body, the lift effectiveness of the wing-body 
ary cidence potentials combination is given by (in Mirels’ notation 
ber Fraenkel’s result for the pressure on an elliptic cone at zero 
ort incidence is in a simpler form than that given by Kahane and I t/a \? a 
Solarski and may therefore be worth quoting. If the equation — a ‘(: ) = (! ) (‘ ) 


of the cone is written 


x = As cos p 
y = Bs sin up 
then 2a/y a2? 
COS T l l 
, 4 1 + (a*/y2? 2 y2? 
Cp = AB| 21n + 
VM? 1(A +B 
1B The same result has been obtained by application of the method 
ae of analysis developed by Adams.? This result does not agree 


A? sin? w + B* cos? wu 
The drag coefficient is 


Cp = AB] 21n : - ] + 0 (A*In? A 
V 


with the equation given by Mirels in an earlier contribution to the 
Readers’ Forum.* However, study of the graphical presentation 
of his results indicates that the latter were based upon Eq. (1 
given above and that the equation given by Mirels in the article 


was apparently incorrectly copied 


od 
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Note on Boundary-Layer Transition in 
Supersonic Flow 


A. H. Lange and R. E. Lee 
U.S. Naval Ordnance Laboratory, White Oak, Silver Spring, Md 
October 15, 1953 


XPERIMENTAL EVIDENCE SHOWS that the laminar boundary 
E layer is made more stable by certain discontinuities in the con 
tour of a body immersed in supersonic flow. This was first noted 
when comparing transition Reynolds Numbers determined on 
cone-cylinders with those obtained on simple cones and on hollow 
cylinders in the NOL 40- by 40-cm. Aeroballistics Wind Tunnel 
No. 2.! The cone-cylinders had consistently shown higher transi 
tion Reynolds Numbers, particularly between the Mach Num- 
bers of 2and3. Transition Reynolds Number is defined by Re = 
vl/v, where / is the wetted length along the model contour be- 
tween forward stagnation point and the point where the boundary 
layer becomes turbulent (determined from spark schlieren pho 
tographs), and where v and » are the free-stream values of the 
air velocity and the kinematic viscosity. The models were 
nearly at equilibrium temperature 

Experiments were recently carried out in the NOL 40- by 40- 
cm. Aeroballistics Wind Tunnel No. 1 with two hollow-cylinder 
models (swallowed front shock) having equal surface roughness 
and either an inside bevel or an outside bevel at the leading edge 
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Fic. 2. Five-degree cone-cylinder model, luminous lacquer coat 
Aeroball.stics Wind Tunnel No. 2, VJ = 3.25 
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Fic. 3 Proposed leading-edge contour of hollow cylinders 


rhe results are shown in Fig. 1 and clearly point out the stabilizing 
effect of the shoulder on the outside bevel hollow cylinder 

The delay in transition may be caused by the favorable pres 
sure gradient that exists for a short length behind a shoulder in 
supersonic flow, although it is surprising to note that the adverse 
pressure gradient, which is weaker but occurs over a larger portion 
of the model, does not entirely override the stabilizing effect 
This stabilization seems to persist even after the boundary layer 
has been completely turbulent for a long distance, as evidenced 
by the luminous lacquer picture of a 5° total angle cone-cylinder 
model (Fig. 2). The picture shows a laminar region on the cy] 
inder (dark) downstream from a turbulent boundary layer which 
extends over approximately two-thirds of the cone 

That a faverable pressure gradient has a pronounced effect in 
supersonic flow can be seen from experiments on circular cylinders 
in transverse flow. In subsonic flow at subcritical Reynolds 
Numbers, laminar separation occurs at a generatrix approxi 
mately 80° from the forward stagnation line, while in supersonic 
flow the separation is delayed to approximately the 115° genera 
trix.” 

If the effect described is of general nature—i.e., not character 
istic for the flow in a wind tunnel only—the question naturally 
arises as to what shape will cause greatest delay in transition and, 
consequentiy, cause minimum skin-friction drag and aerody- 
namic heating. This shape might consist of discontinuities of 
various arbitrary kind, number, and magnitude. In order to ob 
tain some information on this question, tests are planned with 
a family of hollow cylinders having outside bevels of various 
angles but constant frustrum length ( Fig. 3) 
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Transition cn a Standard Model to 


Errata—‘‘Stress Analysis of Sheet-Stringer 


Panels with Cutouts’’! 


Harvey G. McComb, Jr 
Langley Aeronautical Laboratory, NACA, Langley Field, Va 
October 7, 1953 


. HAS BEEN POINTED OUT to the author that certain corrections 


should be made in reference 1 
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1) Each of the coordinates for the bays in Fig. 2(a), page 388, 
should be shifted one bay to the right. Then station (0, 0), for 
tance, will be at the lower left-hand corner of shear panel 
), 0). The corrected coordinates will correspond to those 
hown in Figs. 1, 2(b), and 2(c 
2) The three unknowns in the simultaneous equations writtei 


mn page 396 should be 
Po, Pu, ind g 
instead of 


Px», Pr, and g 





rhe corrected set of unknowns should correspond to the perturba 
tion loads shown in Fig. 6(a), page 397 
3 In the general solution to Eq. (A183), which is the equation 

immediately following Eq Al4), page 399, the first term 
should read 

1 + l 
instead of 

( ry +" l 


ind the second term shou!d read 
A-y—-N y* — 1) 


instead of 


4 The denominator of the bracketed expression in the in- 


tegral in Eq. (A18), p. 400, should read 


3B 2 sin? (x/2Z 


instead of 


3B + 2 sin? (x/2 


This was a typographical error in the manuscript and does not 


invalidate the charts given in the paper 


REFERENCE 
McComb, Harvey G Ir Stre Analyst Sheet-Strineer Panels wit 
tout Journal of the Aeronautical Sciences, Vol. 20, No. 6, pp. 387-401 


June, 1953 


Comment on the Stability Derivatives of a 
Wing-Body-Vertical Tail Configuration 


Arthur E. Bryson, Jr 

Assistant Professor of Mechanical Engineering, Harvard University 
Cambridge, Mass 

ctober 26, 1953 


IT A RECENT NOTE in the Readers’ Forum,’ Mr. Summers 
pointed out that the expression for the apparent additional 
mass of a cylinder with symmetrical wings and unequal vertical 
fins translating along its unsymmetrical axis can be expressed as 
in algebraic function of its dimensions instead of a function con 
taining elliptic integrals as given in my paper.? It can be shown 
that his expression is a simplification of my Eq. (34); my over 
sight in not simplifying the expression was due to the fact that 
the function G*(A, «) defined by Eq. (8) in Appendix (1) can only 
be expressed in terms of elliptic integrals, whereas G(A, uw) = 
G*(X\, wp) + G*(u, \) can be expressed simply as 


G(A, w) = 1 (cos X + cos w}/2 


35) is equal 


implying that the expression in the braces { | in Eq 
to w/2 (which can be shown The expression for the inertia 


coefficient .4;, then becomes (as Summers has shown 
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rhis should have been evident to me from Ward’s paper,* where 
it is shown that the lift involves only the doublet term of the po 
tential, which, in this case, is an algebraic expression, not invol\ 
ing elliptic integrals 

I should also like to correct an error in my paper® in the es 
pression for the inertia coefficient A from which the side fore: 
due to rolling angular velocity is determined, along with severa 
) 


other stability derivatives In Appendix (2) the sign in front of 


the second integral in Eq. (1) should be negative instead of posi 


tive. The correct result for 4 instead of Eq. (36) in th 
paper} is 
2a%h + 
! = sin @ sin # + sin 9 sin @ 
rD 
(h + y2 
sin’ #@ + sin’ @ sin’ @ sin’ @ 
8rD*® (3 
2 cos & cos 4; cos @; cos 6 sin 6 + sin 6 sin @ 
] 
sin @, t cos 6 + cosé cos @ cos 6 


6+6+6 + 6 4 sin 20. + sin 26 
> 


aT 
sin 26, + sin 26, \ 


which is properly antisymmetric if 6) and @; are interchanged with 


6, and @;, respectively For no body (a 0), then @ 8, and 
0; = 4, and A,, simplifies to 
h+f | 
1 = (sin* 6 sin? @ l COS 4 cos 8 
27D (3 
l 
sin @ sin 6, t cos 6 cos 6 
» 
r 1 ‘ 


| ¢ + 6; 4+ sin 26) + sin 24, 
) 
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Three-Dimensional Study of a Jet Penetrating 
a Stream at Right Angles* 


J. P. Fraser 
Engineer, Knolls Atomic Power Laboratory, General Electric 
Company, Schenectady, N.Y 
September 4, 1953 
TWO-DIMENSIONAL TREATMENT, ‘‘Penetration and Defle: 
tion of Jets Oblique to a General Stream" by Fredric | 
Ehrich,! appeared in the February, 1953, issue of the IAS Jour 
NAI A complementary three-dimensional study,? much less 
elaborate, of a round jet penetrating at right angles 1s reviewed 


for comparison 


wresented at Cornell University rhe 
W 


* Based on a master’s thesis 


author gratefully acknowledges the valuable advice and criticism of Dr 


R. Sears and Dr. 5S I. Pai 
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Fic | Distortion of jet cross section in first diameter of travel 

















~| 


x D all 


Fic. 2. Deflection of jet at plane of symmetry 


This study sought to determine the relative importance of 
dynamic pressure forces and mixing to jet penetration. The 
comparison was accomplished indirectly by analytically calcu- 
lating the distortion and deflection of the jet due to dynamic 
pressure forces in the first diameter of travel, extrapolating the 
deflection curve to permit comparison with an experimental 
curve of jet penetration, and inferring that jet penetration is 
chiefly controlled by pressure forces 

The flow within the jet was assumed to be frictionless, irrota 
tional, and incompressible. At the egress plane, jet transverse 
velocity components were assumed negligible. The velocity 
potential within the jet using cylindrical coordinates must 
satisfy the potential equation 

l ] 


Pre T ~('dbr)r TF ooo = UO (1) 
r r? 


The velocity potential within the jet was assumed in the form 






EXPERIMENTAL CURVE 
FIG. 6G, REF 4 





[PRESSURE ESTIMATE 

; > 

4 

Fic. 3. Comparison of calculated and experimental jet de 
flection. The circular path labeled ‘‘pressure estimate’’ is based 
on assuming that one-stream dynamic head acts across an un 
changing jet diameter. Actually, since the pressure difference is 
somewhat greater than one-stream velocity head and the thick 
ness of the jet rapidly diminishes as shown in Fig. 1, the initial 


radius of curvature of the jet corresponding to this analysis would 
be smaller 





o= S + a (7 y" cos n6 + Ux 4 
, i © 
m=1 n=0 
In order to satisfy Eq. (1), the term 4A,,"(7) takes on the form 
Aw? ~ Cr"? 3 
where the C,,” are constants 


The pressure distribution about the jet at the plane of egress 
was assumed identical to that measured about an infinite cylinder 
of the same diameter in a stream of air at the same Reynolds 
Number. This boundary pressure distribution was expanded 


in the series 


Pm Ar 4 


A, cos 86 + As cos 26 4 1, cos 30 4+ 


A, cos 46 +4 j 


> 


with constants determined for streamlines 1, 2, 3, 5, and 7 of Fig 
1. Assuming that cross velocity components in the jet are 
negligible at the egress plane and linearization permissible, pres 


sure may be expressed as 


oe » » ds 
eB 5 
l U' Ox x=0 
By limiting 7 to 4 and m to 2 in Eq. (2), an approximate solu 


tion for the potential function near the egress plane was obtained, 


“, (It is found 


involving five unknown constants, C,", C)}, C 
that, in this approximation, (2° = C,!' = C,4 = 0 These con 
stants were evaluated by Eqs. (4) and (5), at the five stream lines 


mentioned. The result may be denoted by 


Using Eq. (6), radial and tangential curvatures of jet surface 
stream lines were calculated at x = 0. These curvatures were 
used to predict approximate values of 6, 7, ™, M2, and u; at x = 
D/2 for surface stream lines 

The boundary pressure distribution at x = D/2 was based on 
the assumption that the pressure along a surface stream line 
would remain constant for this short distance. This assumption 
appears justified at @ = 0° and in the wake. It allows for the 
spreading of the pressure region on the upstream side qualita 
tively consistent with the flattening of the jet front which takes 
place. 

The flow solution represented by Eq. (6) was abandoned at 
this point in favor of the more elaborate expression 

On 4 = 2(X, r, 0) “> 
m 3 
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satisfying Eqs. (1) and (2) and involving ten constants to be de 
termined 

[he number of unknowns was reduced to five by matching 
static pressures at the same five stream lines as before, at x = O 
Instead of matching the static pressures at more positions at this 
station, the remaining five constants were determined by match 
ing the total pressures at the same five stream lines at x = D/2, 
issuming uniform total head throughout the jet. Actually, to 
perform this calculation, it was found convenient to make use of 
the approximate values of 6, 7, 4), 4 and u; at x = D/2 obtained 
ibove, and to assume the true velocity components to be only 
approximate ones—i.e., to make a 


shortcut, 


slightly different from the 
linear-perturbation calculation. This was simply a 
however, and does not affect the method in principle; the main 
point is that the ten coefficients of Eq. (7) were determined by 
these matching conditions 

Figs. 1, 2, and 3 represent the jet flow corresponding to Eq. (7 
The conclusion was drawn that pressure forces limit jet penetra 
tion by accomplishing the sharp turning after egress followed by 
i further slow penetration due to turbulent mixing. This physi 
cal picture is in keeping with the experimental observations of 
Hawthorne, Rogers, and Zaczek,* that a hole elongated in the 
stream direction gives better penetration than a circular hole of 
the same area; also that penetration is increased if another hole 
or rod is placed immediately downstream 

rhe qualitative agreement between results of this analysis, the 
results of the analysis of Ehrich, and the experimental data re 
ported by Callaghan and Ruggeri may be ¢ xplainable in terms of 
two actions 

(1 The rapid flattening of a round jet as shown in Fig. | mak 
ing a two-dimensional treatment quickly applicable 

2) The filling of the narrow (3 diameters 
by Callaghan and Ruggeri due to flattening and venturi action 
This action is believed ca 


cross section used 


between jet and passage side walls 


pable of contributing toward two-dimensional test data 
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On the Propeller Discontinuity 


L. Talbot and A. K. Oppenheim 
Assistant Professors, University of California at Berkeley 
October 5, 1953 


TT ONE-DIMENSIONAL ACTUATOR DISC CONCEPT has been used 

by several investigators! * for the analysis of compressible 
flow through propellers. The basic assumption common to all 
these studies is that the flow from states 0 to 1 and from states 
2 to 3 is isentropic (the states referred to are indicated in the sys 
tem diagram, Fig. 1, which is drawn for the subsonic flow case ) 
However, various hypotheses have been advanced concerning the 
state change 1-2 through the disc itself. Krzywoblocki' first 
assumed that the specific volume does not change through the 
dise and then extended his analysis to include arbitrary values 


of the specific volume ratio %/7;. Vogeley? and Kitichemann 
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and Weber, on the other hand, considered that the ideal flow 
through the disc is isentropic and that any deviations from this 
process are due to blade friction losses It has alre ady been re 
marked* * that the 


to state 2 is incompatible with the notion of a discontinuous proc 


assumption of isentropic flow from state | 


ess across the actuator disc, and this note is presented in order to 


indicate the consequences of such a discontinuity concept. In 
particular, it becomes clear that entropy increase through the 
disc, which has been ascribed previously to blade friction losses, 
can occur in a frictionless propeller and is imherent in the actuator 
disc concept 

A discontinuity in a flow field is defined as a process that occur 
at constant stream cross-section area and mass flow rate In the 


one-dimensional treatment of the flow system of Fig. 1, it is 


therefore restricted by the following conditions 


Continuity: 


VWomentum: 


Energy: 


defined in Fig. 1; y is the specific heat 


area or dise loading, and ad 


where the symbols are 
ratio, 7 denotes the thrust per unit 
is the power input per unit area of the disc The above equations 
define completely the change of state across the discontinuity, 
provided the relationship between P and 7 is known For the 


incompressible case this relationship is simply 
P=TV,=TV 


However, for the compressible case, the determination of the 
power input involves the integration of the product of local point 
values of the thrust and velocity over the control volume defining 


the discontinuity. For the unidimensional analysis, the result 


of the integration can be expressed as the product of the total 


thrust and some average velocity. An especially interesting 


(1), (2), and (3) results from the simplest as 


solution of Eqs 


sumption for this average velocity—namely, 


V. = (V, + V:)/2 (4) 


The energy equation, Eq. (3), 


, Vite, | - , 


From the momentum equation, Eq. (2), the term in the bracket is 
simply the pressure difference p pi, so that with the use of the 


continuity equation the final result is 
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L (pw2 — piri) = ae = (py pi (6) 
y-1 2 

Eq. (6) represents the locus of possible states 2 for given state 1 
and may be recognized as the Rankine-Hugoniot relations for a 
plane shock wave. The cycle of processes can be represented 
therefore in a form shown in Fig. 2 
the foregoing formulation of the problem leads to a locus of states 
equation which is independent of the thrust and power loadings 


of the disc . oo OE 


Having obtained the locus of states equation for the disc, a 


It is of interest to note that 


specific solution yielding state 2 for a given state 1 can be ob- 
tained by a simple graphical method shown in Fig. 3. The con 
struction follows directly from the momentum equation, Eq 
(2), which, for the purpose of the nondimensional coordinates of 


this figure, can be written as 


(p2/pi) — (1 + 7) 


= 7M," (7) 
1 —- (v2/7) 
where +r = 7/py. 
Since the processes 0-1 and 2-3 are isentropic and the condition 
of equal pressure p; = fo is inherent in the definition of the 


boundaries of the flow system, the change of state 1-2 deter- 
mines explicitly the overall density change 
[(y¥-1) 


P3/ po = € ¥) [(s2—s1)/R] (8) 


where s is the unit entropy and R is the gas constant 
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The flow problem is thus solved with respect to the parameters 
of state 1. However, a propeller problem is usually given in 
terms of the parameter of state 0 (po, %, Vo) and the disc loading 
T. Hence, one more relationship is needed in order to relate 
this state with the discontinuity. This is given by the momen 
tum equation for the control surface 0-3 of Fig. 1—i.e., 


T = mV. V, (9 

A computational solution might proceed as follows: With 
known po, v, and Mo, a value of .W, is assumed and the state 
properties at | are calculated from the isentropic relations. The 


procedure described in Fig. 3 is used to obtain state 2, and the 
cycle is completed with an isentropic expansion back to p; = p 
The computation is repeated 


The value for V; is thus obtained 


with different assumed J1/; until the obtained value of V3 satis 


fies Eq. (9). Calculations for several cases of interest are now 
being carried out by J. D. Stewart and will be presented shortly 
as part of a University of California M.S. thesis 

It is interesting to note that, although all the diagrams have 
a propeller advancing at subsoni 


been drawn for the case of 


speed through still air, the discontinuity analysis presented is 


equally valid for supersonic flow at state 1. Indeed, actual 


computation for a propeller advancing at supersonic speed is 
greatly simplified in that state 1 is the same as state 0, and states 
2 and 38 can be calculated directly without the recourse to the 
iterative procedure described above 

As a final remark, it should be stated that all the arguments 
presented here are restricted by the assumption that the flow 
This definite 


limit on the applicability of the results, since with excessive dis« 


can be treated as one-dimensional imposes a 
loadings the spread of the disturbances in a plane normal to the 
direction of motion and the generation of vortices cannot be 


neglected 
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Calculated Amplified Oscillations in the Plane 
Poiseuille and Blasius Flows* 


S. F. Shen 


Department of Aeronautical Engineering, University of Maryland 
College Park, Md 


October 5, 1953 


OR ACTUAL CALCULATION of hydrodynamic stability problems, 
there are now two somewhat different procedures, one due 


to Heisenberg! and Lin? and the other due to Tollmien’ anc 
Schlichting. 
“inviscid 


The main difference lies in the handling of the so 
called of Blasius 
flow, Lin's neutral curve agreed well with that of Tollmien but 


solutions.”’ In the classical case 
not with Schlichting’s calculation following Tollmien’s procedure 


From the comparison in Schubauer and Skramstad's report,? it 


* This work was done at the suggestion of, and under the supervision of 
Prof. C. C. Lin while the author was at M.I.T. during the summer of 1953 
It was supported by project N5 ori-07872 
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Fic. 4. Comparison of amplification rates (Blasius case). 


appears that the experimental points are somewhat closer to Lin's 
neutral curve than to Schlichting’s. Furthermore, in other de 
tails, such as the amplification rate, Schlichting’s theoretical re 
sult could only be qualitatively verified by the experiment. But 
because of the approximate nature of the theory and the difficulty 
of reproducing the theoretical model in actual experiment, one 
is in no position to draw hasty conclusions on only a single evi 
dence. Recently, Thomas® computed, with the help of the high 
speed IBM machine, several eigen-values of the Sommerfeld-Orr 
equation for the case of plane Poiseuille flow, using finite differ 
ences without resorting to other approximations. Thus we 
now have also a limited number of ‘‘exact’’ results to check against 
the asymptotic methods. Thomas’ values are all complex; it is 
therefore desirable to compute points off the neutral curve for 
direct comparison. The points that lie in the unstable region of 
the (a, R)-plane give, meanwhile, a measure of the amplification 
of the corresponding disturbance. Aside from possible verifica 
tion by available experimental data, the amplification rates also 
might serve as a basis for the development of some criteria to ex 
plain the laminar-turbulent transition, such as the different sugges 
tions of Schlichting,‘ Liepmann,’ and Lees.’ 

We have completed calculations of the amplified oscillations 
for both the plane Poiseuille and Blasius flows by perturbing the 
neutral curve obtained from Lin's procedure. The perturbation 
scheme is the same one previously used by Schlichting in a similar 
calculation. The rates of change of the complex wave velocity 
(« Cr + ic;) in both and a- (wave number) and R- (Reynolds 
Number) directions at the neutral curve are determined. At 
constant R, a cubic in a@ is then fitted as an approximation for the 
intermediate points in the unstable region. Contours of con 
stant c; (representing the amplification rate) are interpolated 
Needless to say, the neutral curve corresponds to c; = 0. The 
results are presented as Figs. 1 and 2. 

For the Poiseuille case, comparison with Thomas’ results are 
made. Remembering the highly compressed scale in the R 
direction (Fig. 1), we note that, except near the critical Reynolds 
Number, the error probably should be estimated at the same value 


of R. It is then seen that a shift in the a-direction of roughly 2 
or 3 per cent is sufficient to account for the discrepancy between 
the values of c;. The critical Reynolds Number by Lin's pro 
cedure, however, is roughly 10 per cent too low Generally 
speaking, the agreement seems to be as good as can be expected 
The upper branch of Lin’s neutral curve might be slightly too 
high (in the a-direction) and the minimum Reynolds Number a 
little too low. Since the method of solution is only asymptoti 
cally correct for large R and since Lin's procedure leaves out 
terms of higher orders in a, the accuracy naturally should be 
somewhat poorer in the lower R- and/or higher a-regions of the 
(a, R)-plane Nevertheless, from the comparison it does lend 
confidence to the adequacy of the asymptotic method 

For the Blasius case, we again make a comparison with Schu 
bauer and Skramstad’s experimental result. To bring out more 
contrast, the neutral curve is replotted in the (8,v/ Uy", R,)-plane, 
where 6,v/ UU? is the dimensionless frequency of the disturbances 
and R, is the Reynolds Number based upon the displacement 
are included in 


thickness. Both Lin’s and Schlichting’s curves 
Fig. 3, together with the experimental points. To compare the 
amplification, the dimensionless amplification rates at R; = | 
630, 1,840, and 2,200 are plotted in Fig. 4* versus the dimension 
less wave number a, again based on the cisplacement thickness 

In both Figs. 3 and 4, the present results are in good agreement 
with experiment, certainly much closer than Schlichting's calcu 
lation 

A detailed account of the work is now under preparation 
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* It might be pointed out that, in reducing their observed data, Schubauer 
and Skramstad used the wave velocity as the downstream propagation 
velocity of the disturbance instead of the more logical ‘‘group velocity.”’ 
rhe correction is to add approximately LO to 20 per cent on the experimental 
amplifications of Fig. 4, thereby bringing a closer check with our computa 
tion at the two higher Reynolds Numbers However, the computed curves 
are obtained by interpolation, which might easily err by as much as 20 per cent 
except in the immediate neighborhood of the neutral curve Hence, we do 


not put too much emphasis on this correction 


Comments on Location of Detached Shocks 
and Truitt’s Dilemma 


E. V. Laitone 
Associate Professor, University of California at Berkeley 
October 5, 1953 


N REFERENCES | AND 2, Truitt mistakenly believes he has de 
rived a new relation for predicting the effect of nose angle on 
the location of a detached shock wave ahead of a plane wedge. 
Furthermore, in reference 1, Truitt claims Moeckel's* result is 
wrong by stating: “ other (i.e., Moeckel’s*) simple expres- 
sions have been derived for shock detachment distance which 


predict no direct effect due to nose shape or varying nose angle 
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rhis prediction ts not true as is borne out by experimental evi 


dence, for thin wedges where the detached shock location is con 
siderably affected by nose angle 

Since Truitt’s result for a wedge is absolutely identical to 
Moeckel's, which does not predict any ¢ ffect of nose shape, 1t 1s 
evident that Duff* was correct in stating: ‘‘(Truitt!) was con 
This is further 
The 


had not been attempted for bodies of revolu 


fused as to the contents of Moeckel’s paper 
borne out by Truitt’s concluding statement in reference 1 
present analysis 
tion at the time of writing.’’ All Truitt had to do was glance 
down two lines in Moeckel’s* paper to find the identical expres 
sion for a body of revolution [Eq. (3) of reference 3 

Now in reference 2 Truitt still claims to have derived an orig 
inal expression ‘which is an explicit function of the opening angle 
of the wedge,”’ showing he is still in a dilemma, although he nox 
admits Moeckel’s 
or affect of nose shape in Moeckel’s* result (or in references 1 or 


result is correct! There is no explicit effect 


2), since the detached shock location is calculated from the body's 
sonic point (as first pointed out by Busemann and Guderley 

g., reference 5) which, of course, is at the shoulder of a wedge 
rruitt’s measurement from the nose of the wedge is merely an 
irtificial shift of the abscissa to a misleading location, as pointed 
out in reference 6, and does not predict any effect of nose 
shape 

One should also apply Duff's! comment to Truitt’s compre 
hension of references 6 and 7, as evinced by his remarks in refer 
ence 2. Reference 7 agrees with Busemann and Guderley that 
the characteristic distance is generally that measured from the 
foremost point of the detached shock wave to the sonic line 
(which must start at the shoulder of a wedge or cone For this 
case the shape of the nose is relatively unimportant (see reference 
5 However, the original claim of reference 7 was that this is 
true only when the sonic line can influence the detached shock in a 
as the Mach 


Number approaches unity, which is the entire basis of the deriva 


finite distance. This is not necessarily the case 


tion presented in reference 7. Therein it is stressed that the de 
velopment provides a limiting value that is restricted to the cas¢ 
where the detached shock and the sonic line become essentially 
two nearly parallel lines, far apart, and extending an unlimited 


distance from the body (see Fig. 2 of reference 7 


This corre 

sponds to Busemann’s® Fig. 9 in the limit as the wedge angle re 

mains small but finite while the Mach Number approaches unity 
so that in the hodograph plane the shock polar reduces to its 
limit at the sonic point. Consequently, it is not surprising that 
a Mach Number of 1.10 may be too large for these limiting con 
ditions to apply. As a matter of fact, all acquainted with the 
transonic problem have noted that the detachment angle pre 

dicted by the transonic similarity relations diverges rapidly from 
the exact value for Mach Numbers greater than 1.10. Obviously, 
the limiting prediction for the detachment distance in references 
6 and 7 would be even further restricted because of the assumed 


independence of the detached shock from the sonic line 


However, the derivation in reference 7 still remains the only 
simple expression having any explicit effect of nose shape, even 
though it rightly may be restricted to Mach Numbers less than 
1.01! For Mach Numbers greater than this (that is, when the 
results of reference 7 do not apply), there cannot be any important 
effect of nose shape, as clearly stated in reference 5. In this cas« 
the only contribution of the nose shape is in locating the sonic 
point; for a wedge this is always at the shoulder. Consequently, 
the usual finite detachment distance is not dependent on the 
wedge angle. 


As a final point it should be mentioned that wind-tunnel tests 
at Mach Numbers approaching unity always place the detached 
shock farther from the body than would be the case in free flight 
because of the shock reflection at the tunnel walls. This type 
of flow would be even more dependent on the sonic line location 


and therefore should follow Moeckel’s’ relatio~ 
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Laminar Boundary-Layer Temperature 
Recovery Factor on a Cone at Nominal Mach 
Number Six* 


Douglas S. Mackay and Henry T. Nagamatsu 
California Institute of Technology, Pasadena, Ca 
October 1, 1953 


A’ INVESTIGATION WAS UNDERTAKEN 1n the spring of 1958 
to determine experimentally the local temperature re 
] 


covery factor for a laminar boundary layer on an insulated 


cone This recovery factor is defined as the ratio of the differ 


ence between the surface and local (just outside the boundary 


layer) temperatures to the difference between the stagnation and 
; I g 


local temperatures. Previous experimental investigations of the 


temperature recovery factors for laminar boundary layers on 


cones in high speed airflow have yielded the results shown in 


Table | 
rABLE | 
Local 
Mach Reynolds Recovery 
Investigators Number Number Factor 
Wimbrow 1.5 3.8 & 10° ) 845 + 0.008 
2 0 9 7-4.8 X 10 oO S55 + 0 OOS 
Eber? 0.88-4.65 6.000-5 * 10 0 845 + 0.008 
des Clers; Stern 
berg 2.18 3.1-7 X 10 ) 851 + 0.007 
Stine; Scherrer* 2.0 +X 106-4 K 10° 0.845 


These results are in close agreement with the theoretical pre 
diction that the recovery factors for a laminar boundary layer are 
essentially constant and independent of Mach Number and Reyn 
olds Number.’ It is the purpose of the present investigation to 
boundary-layer temperature recovery 


thus 


determine the laminar 
factor on a cone for a Mach Number of approximately 6, 
extending the existing data in the higher Mach Number range 

the 5- by 5-in. GALCIT 
continuous-flow, closed-circuit, Hypersonic Wind Tunnel (Leg 
No. 1) with Mach Number 6 fixed-nozzle blocks installed 


20-deg. cone models were used 


The investigation was conducted ir 


Two 
One consisted of a cast ceramik 
core (two parts stabilized zirconia to one part Sauereisen fillet 
and cement), with thermocouples at the surface covered with a 
coating of sprayed steel of 0.010- to 0.015-in. estimated thickness 
(cf. Fig. 1 
thickness with thermocouples soldered to the inside of 


The other was a copper shell of approximately 
0.005-in 
the shell 


eter-Pyrometer, calibrated in degrees Fahrenheit, was used for 


A direct-reading, self-balancing Brown Potentiom 


* The work was sponsored by the Army Ordnance and Air Force 
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THERMOCOUPLES INSIDE 
OF SHELL AT INDICATED 
POSITIONS 


FIG. | 





— THERMOCOUPLE 
LEAD-IN WIRING 


TEMPFRATURE RECOVERY CONE 





The 


temperature was automatically recorded every 45 sec. by a 


the cone surface temperature measurements stagnation 


Honeywell-Brown stagnation con 
Data 


, and pressures, ? 


temper ature 
the 


Minneapolis 


troller-recorder were obtained for flow stagnation 


) 


temperatures, T , and are given in Table 2 


TABLE 2 
P 
Lbs. per Sq.In 
To (°F Absolute ) Remarks 

270 115.2 One phase flow 

267 95.6 One phase flow 

218 95.5 One phase flow 

260) 47.4 One phase flow 

212 17.5 One phase flow 

2138 15.55 One phase flow 
198 95.6 Boundary between one- and 

two-phase flow 

117 95.6 Two phase flow 

108 95.6 Two phase flow 
The free-stream Mach Numbers for the tests are from 5.58 to 
5.88. The corresponding local (just outside the boundary 
layer) Mach Numbers are between 4.78 to 4.96. The free 


stream Reynolds Numbers, computed on the basis of free-stream 
density, velocity, and viscosity, and the lengths along the cone 
ray to the thermocouple locations range between 21,000 and 
540,000. The local temperature recovery factor was determined 
The Prandtl 
Number of the airflow can be computed by the method of refer 


to be 0.844, with an experimental scatter of +0.008 
ence 6. The square roots of the Prandtl Numbers evaluated at 
the average temperatures of the cone surface and of the air at the 
outer edge of the boundary layer were 0.887 and 0.873, respec- 


tively. From the present investigation, the following results 


The local temperature recovery factor for 
0.008 for free 


This laminar recovery 


have been obtained: 
a laminar boundary layer on a cone is (0.844 + 
stream Mach Numbers from 5.6 to 5.9 
factor is independent of the Reynolds Number for the range from 
21,000 to 540,000, and the recovery factor is also independent of 
the Mach Number for Mach Numbers less than approximately 6 
This conclusion is based on a comparison of the present results 
with those of previous investigations at lower Mach Numbers 
The ratio of the temperature recovered on the surface of the cone 
to the stagnation temperature is essentially the same for one 
must be 


phase and two-phase airflows. It kept in mind that 


these conclusions apply only for flow conditions as encoun 


tered in the wind tunnel, where the perfect gas law is valid and 





Where dissociation of air due to high temperatures does not 
occur 
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An Empirical Method for Correction of a Wing 
Downwash Field for the Rolling Up of the 
Trailing Vortex Sheet in Incompressible Flow* 


S. J. Deitchman 


Research Engineer 
Buffalo, N.Y 


October 6, 1953 


Cornell Aeronautical Laboratory, Inc 


pe MECHANISM OF THE ROLLING UP of the trailing vortex 
sheet downstream of a wing, as well as the mathematical 
difficulty of determining the shape of a partially rolled up vortex 
Since for wings of mod- 
a tail 


sheet, has been discussed in reference 1 
erate or low aspect ratio the calculation of downwash at 
location under conditions of partial rolling up might very well 
be required, the following empirical method is proposed for the 
determination of the circulation distribution in the trailing vor 
tex sheet under these conditions. Standard aerodynamic nota 
tion is used throughout, except as stated 

If it is considered that the trailing vortex sheet rolls up into two 
tip vortex cores of finite diameter and constant angular velocity 
throughout,! then on the basis of the Helmholz vortex theorems 
it may be assumed that the tip vortex cores contain all the fluid 
particles carrying vorticity and that these particles were originally 
distributed along the wing trailing edge. A mechanical model 
can be visualized in which all the vorticity-carrying fluid that 
leaves the wing trailing edge at a particular instant is carried in a 
perpendicular to the horizontal axis of symmetry, of 
finite unit thickness in the x direction. At the time ¢ = (0), the 


fluid is considered to be contained in this “plane” in a rectangular 


“plane” 


‘wake”’ at the trailing edge. Since the distribution of fluid par- 
ticles carrying vorticity is arbitrary, the spanwise variation of 
strength of the trailing vortex sheet is accounted for by having 
flow spanwise into the cores, as the 


the fluid in the ‘‘wake” 


“plane” is carried downstream at the air-stream velocity, V, 
with a velocity that varies with the distance from the core center 
Then the problem may be considered one of fluid flow in a channel 
of unit depth and height at each spanwise station proportional to 
the velocity there, starting at time ¢ = (0 

The vorticity in the spiral coils of the trailing vortex sheet 
near the cores is assumed to be carried in the cores. Then the 
initial circulation distribution on the wing is the same as the vor 
ticity distribution in the channel at ¢ = 0, and at any subsequent 
time the effective circulation distribution is given by the dis- 
tribution of fluid between the ‘‘wake"’ and the tip-vortex cores 
The span, 6’, of the trailing vortex sheet at any instant (i.e., the 
distance between centers of the tip vortex cores) may be deter 
mined from the circulation distribution curves by a graphical 
integration, since from the requirement of conservation of im- 
pulse the areas under all the distribution curves must be equal 
The vorticity concentrated in the cores is given by the value, 
r,, of the circulation at the point y = 6'/2. When the rolling 
up process has been completed, the distribution will be rectangu 
lar, of span 6,’, and the circulation will be equal to that in the 
wing plane of symmetry (disregarding small areas of negative 
vorticity there). The process is illustrated in Fig. |! 

In order to map the rolling-up process, it 1s necessary to assume 
a small initial value [,,, at an initial 6’/2, as shown in Fig. 1, and 


a small initial radius, r,., say 0.15 7 The final core radius is 


given by 


$7ACy 
=. exp Sgr 
4/2 h | cA 
as determined from energy considerations,” and the initial wake 
height is determined from r,. and the wing span, 6: 
* Taken in part from the writer’s M.S. thesis written for the University of 
Buffalo 


any point, x = 
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Fic. 1. Steps in the rolling-up process 
2rr = D 2 


For most wing loadings, the spanwise velocity at any pomt, Vo, 1 


the ‘“‘wake"’ may be taken as 


where the velocity entering the core is taken as 


The initial height of the ‘‘wake”’ in the plane of symmetry at / = 


(O)* is determined by spanwise integration as 7.6 2, and 
° - Zs 
the solution of the channel flow problem for incompressible flow 


yields the following results: 


Considering the velocity of a particle constant over a short 
period of time and the ‘“‘wake”’ height at the tip small compared 


with the core radius, the fluid at yo will move in the short time 4, = 


x,/V to 


The core radius will be given by 


o C, ( ) b : | D 
4.0 - r. 2% _ 4. oe 
, AL\c,c/,, 2 V2 


r,%* = 0 . 
6 


Then the rolling-up process may be solved to completion or to 
a)» 
— ‘ 


, by numerical integration over a series of 
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small steps, «;, wherem the results y,, 0,’ and r,, of each step are 


considered the initial conditions, r,,, bo’, and y, of the succeed 
ing step. The results have been given in terms of the commonly 
used parameter 

ic/C,é = 2AP/bVC, 


l 


Cc 


The reasoning of Westwater® concerning the assumption of two- 
dimensionality of the trailing vortex sheet at any instant may be 
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used to justify downwash calculations at the tail using the circu 
lation distribution determined for the partially rolled up trailing 
vortex sheet at the tail location. The method should be appli 
cable to wings with sweptback, as well as straight, trailing edges 
if it is considered that no rolling up occurs until the trailing vor 
tex sheet leaves the wing tips. Many experimenters! * * have 
found that the tip vortex cores will move downstream at about 
the level of the wing tips, while in the plane ot symmetry the 


wake displacement will be given by 


ex 
h = { tan € dx 
Te 


. 


(see reference 5 There exists in the literature sufficient water 
tank data for a rough estimation of the wake displacement at 
other spanwise locations 

Calculations have been made of the rolling-up process for an 
elliptical wing of A = 5at a = 8°, the results of which are shown 
For comparison, the locations of the tip vortex cores 
Since 


the square root law assumed for the spanwise velocity distribution 


in Fig. 2a 
determined from water tank experiments! are also shown 


in the ‘“‘wake”’ is most valid for elliptical or near-elliptical circu 
lation distributions, an extreme case was tested using the roughly 
triangular distribution given for a delta wing of A = 2.04 ata = 


» 


33° in reference 1. The results, shown in Fig. 2b and compared 
again with water tank data,! show that, even in this case, the pres 
ent method will give an answer that is qualitatively acceptable for 
design purposes 

For a rough estimate of the distance to complete rolling up, the 
equation? 
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eful \ ilues ol A have been determined for circulation dis 
tributions from rectangular (A = 0) to triangular and are given 
in Mig. 3 
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Some Properties of Intense Sound Waves in 
Relation to Combustion 


J. C. Truman 
Assistant Professor, Aeronautical Engineering Department, USAF 
Institute of Technology, Wright-Patterson Air Force Base, Ohio 


October 26, 1953 


A* IMPORTANT PROBLEM IN COMBUSTION is an understanding 
~ of the interrelationships among unsteady flow conditions 
and various combustion phenomena. This problem becomes 
particularly significant when one considers a combustion system 
in which the flow is rendered inherently unsteady by the design 
of the system or in which combustion instability or rough burning 
may occur. An approach to the problem has been to investigate 
experimentally the effects of sound waves on various combus 
tion parameters; in several instances, sound waves of relatively 
high intensity have been used. It is the purpose here to sum 
marize from the literature certain second-order effects that are 
characteristic of intense sound waves (as differentiated from 
single shock waves on the one hand, and periodic but small acous 
Some of these effects may be signin- 


laboratory 


tic waves on the other 


cant in the interpretation of data from tests on 


simple burners and in the analysis of combustion instability in 
more complicated systems 

The distortion of the wave shape of intense sound waves in air 
has long been recognized. It was treated by Rayleigh! and 
Lamb,? who also referred to the earlier works of Poisson, Stokes, 
Earnshaw, Riemann, Rankine. Briefly, this distortion 
stems from the nonlinearity of the equations of motion. Con 


sider the pressure-displacement plot of a finite amplitude sinu 


and 


soidal wave. Since the high-pressure portions of the wave travel 
faster than the parts having lower pressure, the effect is a steep 
ening of the wave front as it propagates. Since physical con 
siderations require that pressure be single-valued at any point, 
i limit exist; the limiting 


mechanism is provided by viscosity 


to this steepening process must 
An analysis by Fay,* for the case of plane waves, shows that 
nonlinearity results in the generation of harmonic components 
as the wave propagates. Also, extraneous components having 
sum and difference frequencies are established by the interaction 
of harmonics. The result is a continuous transfer of energy from 
components of lower frequency to those of higher frequency 
The effect of viscosity is to attenuate the higher frequency com 
ponents more than the lower. Thus, an equilibrium wave form 
exists, in which viscous attenuation of higher frequency com 


ponents counteracts the energy transfer to those components 
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All waves approach this equilibrium form, which depends on fre 
queney and amplitude It is not permanent, but changes from a 
sawtooth for a wave of large amplitude to a sinusoidal shape as 
amplitude becomes infinitesimal 

Phuras, Jenkins, and O’ Neil‘ have investigated the generation 


of harmonics, deriving expressions for the pressure amplitudes 


of the second and third harmonies as functions of the distance 
from the source and the frequency and pressure amplitude of the 
highet 


Other 


fundamental. As these parameters are increased, the 
harmonics are found to become relatively more important 
contributions to the theoretical study of the problem appear in 
references 5-11 


Experimental evidence of wave distortion is given in refer- 


Actually, the data presented in reference 16 were 


ences 12-16 
obtained in a rocket motor during combustion, and the observed 
distortion is not necessarily attributed entirely to nonlinear ef 
fects 


Other second-order effects that have been observed in intense 


sound fields are an excess static pressure!* and the so 
called ‘‘sonic wind,’’ a d.c. flow of air away from the sound 
source. !4) 14 18 These effects are relatively small in magni 


tude and presumably are not important in most combustion work 


But distortion and the generation of harmonic components 


may be highly significant in the investigation of combustion in 
stability in rocket motors and other similar types of burners 
These effects may explain some of the frequencies that have been 


measured in such investigations 
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On Alternative Forms for the Basic Equations 
of Transonic Flow Theory 


John R. Spreiter 

Aeronautical Research Scientist, Ames Aeronautical Laboratory 
NACA, Moffett Field, Calif 

October 29, 1953 


A’ TENTION HAS BEEN CALLED by numerous authors! ? to the 
possibility of certain alternative forms for the equations for 
transonic flow about thin wings. It is the purpose of this note 
to contribute to this discussion and to indicate some reasons for 
the selection of one form of these in preference to another more 
widely used form 

Consider a rectangular coordinate system xyz so oriented that 
the x axis is parallel to the free-stream velocity; let ¢ represent 
the velocity potential and ¢ the perturbation velocity potential 
so that 


y= Ux (1) 
If it is assumed that the local velocities differ only slightly from 
the free-stream velocity, an approximate differential equation for 


y appropriate for the transonic speed range is 
(1 Mo )er2 + Oyy + ¢22 = Moe? Or¢rz = kereze (2) 


Woy + 1)/Uo. The corresponding approximate 


relation for the shock polar is 


where k = 


(? = MPV lee, — en}? + hen — On? 


/ ° ” 
(Px, ro)” (od) 


where the subscripts 1 and 2 represent the conditions on the two 
sides of the shock wave. Upon solving the above equations sub 
ject to the appropriate boundary conditions for thin wings, we 
may calculate the pressure coefficient to the same approximation 


using the following relation: 

Cp = —(2/Uo)exz (4) 
It should be noted that the results obtained by using the above 
equations might be expected to tend toward those of linear theory 
as the free-stream Mach Number MM, departs far from unity. 
This follows from the fact that the product g:¢,, becomes small 
relative to the linear terms under this condition. Solutions of 
the equations for transonic flow found to date have all possessed 
this property. 

In common with the pioneer work of von Karman,’ most of the 
transonic results published in the U.S.A. have been based on a 
different set of equations obtained by assuming that the local 
velocities differ only slightly from the critical speed of sound a,. 
Thus, select the coordinate system the same as above; let ¢ again 
represent the velocity potential and g’ the perturbation velocity 
potential so that 


¢ = @ — aX (9) 


Note that, in this formulation of the problem, the perturbation 
velocities do not vanish far upstream when the free-stream Mach 
different 
generally used for these quantities are 


Number is from unity. The approximate relations 


, 9 a, , , 
(gr’)o = (Mo? — 1) E (gy lo = (¢:')o = O (6) 
: ae 1 
The approximate differential equation for ¢’ is 


; Pm or: ak 


, , ~ 
Puy G2 = Er Orzx (a) 


ay 
The corresponding relations for the shock polar and the pressure 


coefficient are 
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In the application of results obtained from the latter set of 
equations, it has been noted that the theoretical results appear 
Mach 


Number is far removed from unity, in spite of the fact that the 


to be applicable quantitatively even if the free-stream 
equations from which they are derived are valid only in the im 
mediate neighborhood of sonic speed. The reason for this be 
came clear when it was demonstrated‘ that the values obtained for 
C,, and hence for integrated forces and moments, using Eqs. (5 
through (9) are identical with those obtained using the less re 
strictive Eqs. (1) through (4), provided Wy is equated to unity in 
2) and (3). Thus the coeffi 
cient & is approximated by (y + 1)/ Uo rather than WW o%(y + 1 

ly. This might seem like a reasonable approximation, inasmuch 


the right-hand members of Eqs 


as these particular terms are merely approximations to allow the 
treatment of transonic flows and rapidly diminish in magnitude as 
JJ, departs from unity. It will be shown in the remainder of this 


note, however, that the retention of the full coefficient as given 
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in Eqs. (2) and (3) is of considerable importance and leads to a 
superior formulation not only from the point of view of theoretical 
considerations but also from that of actual comparison of theo 
retical and experimental results 

Consider first the prediction of the variation with free-stream 
Mach Number of the critical pressure coefficient Cp,, defined as 
the value of the pressure coefficient C, at a point where the local 
Mach Number isunity. Inthe transonic theory, Cp-, corresponds 
to that value of C, and, hence, of ¢,, at which Eq. (2) changes 
from elliptic to hyperbolic type. This condition is recognized 
by the vanishing of the coefficient of ¢,., thus 


| M? 


(10) 
or, in view of Eq. (4), 


(11) 


(Gx Jer 


Us ~~ BUs 


Con 


The exact relation for isentropic flow is (e.g., reference 5, page 28) 


9 2 . 1 7 1) 
Cper > 2 a : : Me) (12) 


The variation of Cp-, with 7) has been computed using the exact 


and approximate relations, and the results are presented in 


Fig. 1. The superiority of the result obtained when 1/,? is re 
tained in the coefficient & is readily apparent 

A second case where the exact and approximate relations can 
be compared is furnished by considering the velocity jump 
through a shock wave. If the flow ahead of the shock wave is 
uniform and parallel to the x axis, the results may be conveniently 
represented by the shock polar diagram in which V2 is plotted as 
The exact relation (e.g., reference 5, page 53) is 


U,U-2 a 
+ 1))/U,? 


a function of ls. 


UUs +a 


* 


where Ul’ and V refer to velocity components parallel and per 


pendicular, respectively, to the x axis. The corresponding ap 
proximate relation is determined from Eq. (3) by setting ¢.,, 


¥y» x, and ¢,, to zero, whereby 


cal Gre” (14) 


The corresponding variation of V; with lL’ may be readily deter 


mined, since, for this case, 


Us. = Uo : y= . (15) 


The variation of 1» with l% for WW) = 1.2 has been computed 
using both the exact and approximate relations, and the results 


are presented in Fig. 2. It may be seen that the retention of 
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7)? in the expression for k provides the exact relation for the 
velocity jump through normal shocks and greatly improves the 
quality of the approximation for all intermediate cases 

In problems such as we are considering here, the final test is 
Although 


both experimental and theoretical results for the transonic speed 


provided by comparison with experimental results 


range are limited, complete information does exist at the present 
time for the drag of a single-wedge section followed by a straight 
section extending far downstream.'» * ‘ These theoretical re- 
sults were determined originally using Eqs. (5) through (9) and 
are therefore identical with those that would be obtained using 
Eqs. (1) through (4), provided k is equated to(4 + 1)/U>. How- 
ever, they may be easily converted to those that would be ob- 
tained using the alternative expression for k by simply replacing 
(vy + 1)/Uo with Wy2(y + 1)/U> wherever it occurs. 

Fig. 3 shows the theoretical and experimental results plotted 
in the same manner as in the original papers. The small vertical 
lines on the experimental data points represent the uncertainty 
of the values. This figure indicates that the theoretical and ex- 
perimental results are only in general qualitative agreement when 
k is approximated by (y + 1)/l 

The same results are replotted in Fig. 4, with k equated to 
Woy +1 It can be seen that the 


theoretical and experimental results are now in nearly perfect 


Uy rather than (y + 1)/l 


agreement. Comparison of Figs. 3 and 4 provides striking evi- 
dence supporting the contention that k should be equated to 
Wo2(vy + 1)/1 than (7 + 1)/l 
matters are discussed at greater length in a forthcoming NACA 
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rather These and related 


papet 
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